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Abstract. The different facets of the i?-matrix method are presented pedagogically 
in a general framework. Two variants have been developed over the years: (i) The 
"calculable" i?-matrix method is a calculational tool to derive scattering properties 
from the Schrodingcr equation in a large variety of physical problems. It was developed 
rather independently in atomic and nuclear physics with too little mutual influence, (ii) 
The "phcnomcnological" i?-matrix method is a technique to parametrize various types 
of cross sections. It was mainly (or uniquely) used in nuclear physics. Both directions 
are explained by starting from the simple problem of scattering by a potential. They 
are illustrated by simple examples in nuclear and atomic physics. In addition to 
clastic scattering, the i?-matrix formalism is applied to transfer and radiative-capture 
reactions. We also present more recent and more ambitious applications of the theory 
in nuclear physics. 
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1. Introduction 

1.1. Principle 

The i?-matrix theory is a powerful tool of quantum physics, introduced by Wigner 
and Eisenbud [U [21 E] where they simplified an original idea of Kapur and Peierls [I] . 
The advantage of their simplification is that the R matrix only involves real energy- 
independent parameters. Initially the theory was aimed at describing resonances in 
nuclear reactions. However even the very first developments also contained the principle 
of a technique for solving coupled- channel Schrodinger equations in the continuum. 

At present, the main aim of the R- matrix theory is to describe scattering states 
resulting from the interaction of particles or systems of particles, which can be nucleons, 
nuclei, electrons, atoms, molecules. Its principle relies on a division of the configuration 
space into two regions: the internal and external regions. The boundary between 
these regions is defined by a parameter known as the channel radius. This radius is 
chosen large enough so that, in the external region, the different parts of the studied 
system interact only through known long-range forces and antisymmetrization effects 
can be neglected. The scattering wave function is approximated there by its asymptotic 
expression which is known except for some coefficients related to the scattering matrix. 
In the internal region, the system is considered as confined. Its eigenstates thus form a 
discrete basis which can be calculated. A scattering wave function at an arbitrary energy 
is expanded in the internal region over these square-integrable eigenstates. Then, the 
R matrix, which is the inverse of the logarithmic derivative of the wave function at 
the boundary, can be calculated. A matching with the solution in the external region 
provides the scattering matrix. This method can also provide the bound states of the 
system. In this case, the external solution behaves as a decreasing exponential. Since 
the exponential decrease depends on the unknown binding energy, an iteration is then 
necessary. 

The i?-matrix theory was developed into two different directions with little exchange 
between these variants. Many of its practitioners often ignore the progresses about the 
other aspect of this double-faced method. 

As already mentioned, the original goal was to provide an efficient theory for the 
treatment of nuclear resonances [31 [5]. From information on bound states and low- 
energy resonances, it soon became clear that the i?-matrix theory offers an efficient way 
for accurately parametrizing not only resonances but also the non-resonant part of low- 
energy cross sections with a small number of parameters [5] . An important advantage is 
that most of these parameters have a physical meaning. This first variant of the method 
is still very important and much employed, in particular to parametrize the low-energy 
cross sections relevant in nuclear astrophysics. This version of the i?-matrix theory will 
be called hereafter the phenomenological R matrix. Its properties are reviewed in [51 Ej . 

The other aspect of the i?-matrix theory is that it can provide a simple and elegant 
way for solving the Schrodinger equation. It is especially competitive in coupled-channel 
problems with large numbers of open channels where direct numerical integrations may 
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become unstable. An additional advantage is that narrow resonances which can escape 
a purely numerical treatment are easily studied. This other facet of the .R-matrix theory 
has been mostly developed in atomic physics although we shall see that it can also be 
very useful for nuclear-physics applications. This variant will be called hereafter the 
calculable or computational R matrix. Its properties are reviewed in [3, [HI EJ ITO] . 

A very comprehensive review of the phenomenological .R-matrix method has been 
given in 1958 by Lane and Thomas [5]. Their article contains most of the important 
aspects of the phenomenological applications of the R matrix to nuclear physics. Many 
of their results can also be useful for the calculable R matrix. However, in 1957, just 
before that review appeared in print, an important improvement of the method was 
published which is therefore not used in their review. Bloch introduced a singular 
operator defined on the boundary between the two regions, now known as the Bloch 
operator, which allows a more elegant and compact presentation of the method [TT] . 
The main interest of the Bloch operator is that its use led to extensions of the method 
to more general treatments of the resolution of the Schrodinger equation in the internal 
region and opened the way to accurate methods of resolution in atomic and nuclear 
physics. Several reviews on the computational R matrix have been published in the 
context of nuclear physics [7] and of atomic physics jHJ [121 HQ] - Reference [9] deals with 
both aspects. An update is nevertheless timely. 

1.2. The phenomenological R matrix 

The phenomenological R matrix and most of its applications were already exhaustively 
described fifty years ago in [5]. Among these applications, let us mention a detailed 
study of resonances and an extension of the method to the description of electromagnetic 
processes. As far as we know, all these applications have been made in nuclear physics, 
i.e. for the scattering of neutrons on nuclei or of nuclei on nuclei with the presence of 
a repulsive Coulomb barrier. Nevertheless, using the method still revealed a number of 
difficulties. In a series of papers, Barker and collaborators provided practical solutions 
to the determination of the .R-matrix parameters from experimental data [131 H31 E] 
and applied this framework to the spectroscopy and reactions of light nuclei [T51 [T7] . 
They also explained non-intuitive effects such as the Thomas- Ehrman shift p2], ghosts 
of resonances [18] and extended the method to further processes such as radiative- 
capture reactions [191 EDI EI] and delayed (3 decay [22J. The approach developed by 
Barker and collaborators has become a standard tool for the analysis of low-energy 
radiative-capture reactions useful in astrophysics. Recent progresses in the adjustment 
of .R-matrix parameters have been performed in [231 EI] • 

The R matrix allows parametrizing various physical processes and its determination 
provides collision matrices and cross sections. For each set of good quantum numbers, i.e. 
total angular momentum and parity, the dimension of the phenomenological R matrix 
is equal to the number of channels relevant to the physical properties. When a single 
channel is considered, the R matrix for a partial wave with orbital momentum I and 
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total angular momentum J is a function of the energy E parametrized by the formula 



In principle, this function possesses an infinity of poles at the real energies E n \ j but only a 
limited number N of such poles affect the low-energy cross sections. The lowest poles are 
closely related to bound states at negative energies or to narrow resonances at positive 
energies. Nevertheless, the poles and the energies of physical states are slightly different. 
Because of this shift, the determination of these parameters from data requires some 
skill. The real parameters are known as the reduced width amplitudes because their 
square is a crucial factor of the width of non-overlapping resonances. More precisely, we 
shall see in section [3761 that the width is given as T = 2j%ijPi where P L is the penetration 
factor which includes most of the effects of transmission through the Coulomb barrier. 
This factor depends on energy and the width thus also depends on energy. 

A serious drawback of the phenomenological R matrix is that the poles and widths 
depend on the choice of channel radius, i.e. on a rather arbitrary value. This aspect 
of the R matrix has been criticized by a number of authors and has led to further 
developments of the competing phenomenological K matrix [25] . The K matrix, which 
provides an alternative formulation of the collision matrix, is also expanded in a series 
involving an infinity of poles. This approach is based on a delicate treatment of Coulomb 
functions. In spite of the fact that the K matrix does not contain an arbitrary parameter 
such as the channel radius, its parametrization is more difficult because its parameters 
may have a less direct physical interpretation. 

1.3. The calculable R matrix 

The aim of the calculable R matrix is to provide an efficient way of solving the 
Schrodinger equation both at positive and negative energies. It was proposed in 1965 
by Haglund and Robson and applied to a two-channel problem involving square-well 
potentials [26]. An expansion over a finite basis was introduced by Buttle [27]. He 
performed the first realistic application on 12 C + n scattering [27]. He also proposed 
a correction to the truncation of the R matrix to a finite number of poles, that is now 
named after him. A more serious problem is a discontinuity of the derivative of the 
wave function at the boundary between the regions that occurs with the traditional 
choice of basis states inspired by the original ideas in [3]. Various solutions to the lack 
of matching at the boundary have been suggested (see [9] for a review). This apparent 
problem has attracted a lot of attention even long after an efficient technique where it 
does not occur was introduced [28j [29] . By dropping an unnecessary condition as we 
will describe, the .R-matrix method can be very accurate without matching problems 
and without need for a Buttle correction. 

Some users of the phenomenological R matrix consider the channel radius as a 
parameter which must be optimized when fitting the data. Even if this dependence on a 
parameter without strong physical meaning is weak, this is a drawback that would not 
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be acceptable when aiming at accurately solve the Schrodinger equation. Hence a crucial 
test of the results of the calculable R matrix is an almost perfect independence with 
respect to the choice of channel radius. This test provides a measure of the accuracy of 
the calculations. 

In spite of its introduction for nuclear-physics problems [251 121], this approach 
was first extensively developed to study electron (or positron) collisions on atoms and 
molecules [El 121 EI]- It allows describing the excitation and ionization of these systems. 
Photoionization, i.e. collisions with a photon leading to the single or double ionization 
of the atom, is also a well-studied application [30| I3T]. 

Important and difficult aspects of these atomic-physics problems are the non- 
locality of the interaction due to electron exchanges and the long-range nature of the 
interactions, due to the polarization interactions. The non-locality is well treated in the 
.R-matrix approach. The long range of the force implies that the asymptotic behaviour 
of the solution is only reached for very large values of the interparticle distance. To 
avoid using a very large channel radius, propagation methods have been introduced 
[32J. They involve an intermediate region where the interaction can be simplified, for 
example with an asymptotic expansion. 

Electron scattering on heavy atoms has required the introduction of relativistic 
corrections. The extension of the R matrix to the Dirac equation has been introduced 
as early as in 1948 [33] but its validity remained controversial during a long time. 
Relativistic corrections were thus first derived from the Breit-Pauli equation. The 
validity of the Dirac extension is now well established and elaborate relativistic codes 
have been developed. This aspect will not be covered here (see [M] for a recent review). 

In nuclear physics, the computational R matrix is much less used although it 
should be very useful in large coupled- channel calculations. It has been much applied 
in microscopic cluster calculations in which the difficult antisymmetrization is taken 
into account in the internal region only [351 EE] . The versatility of the R matrix found 
interesting applications in processes where bound and scattering states are mixed, such 
as radiative capture or delayed j3 decay [371 EH]- The application of the .R-matrix 
method to coupled-channel calculations has been simplified by its combination with 
the Lagrange-mesh method which avoids calculating matrix elements of the potentials 
[391 HOI E]. Recently this approach has been extended to non-local interactions [32]. 

Other approaches to the same problem present a number of similarities. The 
variational A'- matrix method [13] has a very similar spirit. The Kohn variational 
principle for the logarithmic derivative is equivalent to the calculable R matrix [H]. 
Practical implementations of the Gamow-state method [IS] are also exactly equivalent 
to .R-matrix calculations [46j . 

1.4- Outline 

In this review, we present both calculable and phenomenological versions of the .R-matrix 
theory. Very few papers deal with both aspects simultaneously [UJ HH]. Since many 
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excellent reviews already exist, we try to make an introductory presentation, illustrated 
with simple numerical examples. We also show the parallel evolutions of the method in 
atomic and nuclear physics and try to shed light on some common misunderstandings or 
controversies about the i?-matrix methods. The numerical examples that we display are 
tailored to allow a motivated reader to test his/her understanding by reproducing them 
with limited effort. Because of our background as nuclear physicists, most examples 
(but not all) correspond to nuclear applications. Finally, we also review state-of-the-art 
calculations in nuclear physics where the R matrix proves useful. 

Contrary to tradition, we start with the calculable R matrix on a finite basis, 
which provides a convenient numerical approach. Taking the limit for an infinite 
complete basis will introduce the theoretical R matrix which leads after truncation 
to the phenomenological approximation (11.11) . For the sake of simplicity, we detail 
potential scattering in the single- channel case so avoiding the unpedagogical definitions 
of channel wave functions. We only provide the main steps for the multichannel case. 

The bibliography about the R matrix is enormous and can not be fully covered 
here. We have tried to quote papers that we think significant or useful for further 
bibliographic research. 

In section [2j we introduce the necessary basics of scattering theory with radiative 
capture as a more elaborate application. The calculable R matrix is presented in section 
[3] and relatively simple numerical applications in section HI The phenomenological R 
matrix and and its applications are presented in section [51 Recent elaborate calculations 
in nuclear physics are reviewed in section Concluding remarks are made in section [7J 
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2. Summary of scattering theory 

2.1. Coulomb scattering 

Consider the collision of two particles with respective masses m\ and m 2 and charges 
Z\e and Z" 2 e at a positive energy E in the centre-of-mass frame. The wavenumber is 
defined as 



k = yj2(iE/h, (2.1) 

where /x = m^/ (mi + 7722) is the reduced mass. 

Let us start with some definitions about pure Coulomb scattering. In this 
Bohr radius can be defined as 

a B = ,„ „ , 2 - 2.2 

A useful parameter is the dimensionless Sommerfeld parameter 

= Z~iZ 2 e 2 = sgn(ZiZ 2 ) 3 
/if a#/c 

where u = hk/n is the relative velocity. Parameter 77 measures the importance of 
Coulomb effects at a given energy. The neutral case is recovered with 77 = 0. 

For a central potential, a wave function can be factorized in spherical coordinates 
r = (r, Q) as ip(r) = r~ 1 ui(r)Y l m (Q). The spherical harmonics YJ m (f2) depend on the 
orbital and magnetic quantum numbers I and m, and on the angles Q = (9,(p). They 
are defined according to the convention of Condon and Shortley. The radial Schrodinger 
equation for the Coulomb problem in partial wave I then reads 

/ d 2 1(1 + 1) 2kr] , 9 \ , , , , 

(S5 " " -T + « 2 ) Mr) = 0. (2.4) 

Its solutions are combinations of the regular and irregular Coulomb functions Fi(rj, kr) 
and Gi(rj,kr) [19]. The regular function vanishes at the origin and is normalized at 
infinity according to 

Fi(r),x) — > sm(x — \lix — ri In 2x + aA, (2.5) 

where appears the Coulomb phase shift 

<7 { = argr(Z + 1 + 177) (2-6) 

involving the Euler function T. The irregular function Gi(rj,x) is unbound at the origin 
(except for 77 = I = 0) and is fixed by its asymptotic behaviour 

Gi(r),x) — > cos(x — \lit — 77 In 2x + aA. (2.7) 
Also very useful are the conjugate functions 

7,(77, x) = G l - %F h 0,(77, x) = Gi + iF h (2.8) 
which behave asymptotically like incoming and outgoing waves, respectively, 

7,(77, x) — > e -i(*-b«-v^+<n) ^ 0,(77, x) — > e^-f^-^^). (2.9) 
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Notice that we do not follow here the same phase convention as Lane and Thomas [5] . In 
the neutral case, one has Fi(0, x) = xji(x) and Gi(0,x) = xrii(x) where ji and rii = —yi 
are spherical Bessel functions [49] . 

In some applications, solutions of (12. 4p are also needed at negative energies. The 
real solution decreasing at infinity is the Whittaker function W_ VB l+ i (2nr) |49j . It 
depends on the wave number k = ^—2fiE/h and on the Sommerfeld parameter 
t]b = sgn(ZiZ 2 )/aBK of the bound state. Whittaker functions behave asymptotically as 

W -W+|(*)^-^ /2 . (2.10) 

They are singular at the origin. 

A bounded solution of the three-dimensional Schrodinger equation at a positive 
energy with the same Coulomb potential is given by [50] 

V£(r) = (27r)- 3 / 2 e^ 7r,?/2 r(l + i V )e ikz x F x {-ir), 1, ik(r - z)) (2.11) 

where \F\ is the confluent hypergeometric function [49J. This wave function has the 
asymptotic behaviour of an outgoing scattering state 

e i[k z+V lnk(r-z)] + jtQ, e ^ 

r J 

The coefficient of the second term is the Coulomb scattering amplitude, 

The square of the modulus of fc provides the Rutherford cross section. Function ip^ 
can be expanded in partial waves as 

oo 

^+(r) = (2ny 3 / 2 (kry 1 Y / (^ + ^i l e m Pi(cos 9)^(7],^), (2.14) 

i=o 

where Pi is a Legendre polynomial [49]. 
2.2. Scattering by a potential 

Consider a potential V tending to the Coulomb potential faster than r~ 2 , 

nr)^*^ + o(±). (2.15) 
The radial Schrodinger equation in partial wave I reads 

(£-^-^ + *)«.M = o ( 2 ,a) 

with the condition at the origin 

u,(0) = 0. (2.17) 
A real solution at positive energy E behaves asymptotically as 

ui(r) — ► cos SiFi(rj, kr) + sin SiGi(r],kr), (2.18) 
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up to a normalization factor. The important physical quantity is the phase shift 5i. For 
later use, it is however more convenient to write the solution as 



u i W ^ tfiftfa fcr ) - Wifa, (2-19) 
where C\ can be chosen in various ways. For example, (12.1 8p is recovered with 
Ci = iexp(—i5i)/2 and the normalization of the ui with respect to 5(k — k') is obtained 
with Ci = iexp(—i5i)/^/2ii. The collision or scattering 'matrix' Ui is given by 

Ui = e 2iSl . (2.20) 

With the different partial solutions, one can construct an outgoing stationary 
solution 

oo 

m + (r) = (27ry 3/2 (2kr)- 1 £(2Z + l)i l+l e tai P^cos 6)Cf l Ul (r) (2.21) 

1=0 

behaving asymptotically as a Coulomb wave (12. lip propagating in the z direction plus 
an outgoing spherical wave 

p i(kr—r) In 2kr) 

¥+(r) ^ ^(r) + (2 7 r)- 3 / 2 /(^) (2.22) 

The coefficient of the second term in this asymptotic expression determines the 
additional scattering amplitude 

1 oo 

= ^7 T,( 21 + l)e 2 ^ - l)Pi(cos9). (2.23) 
From (I2.23p . one obtains the elastic cross section 

^ = l/c(fi) + /(fi)| 2 (2.24) 

also involving the Coulomb amplitude (12.131) . The scattering wave function (I2.2ip is 
useful in various types of reactions. We illustrate it below with radiative capture. 



2.3. Collisions in a many-body system 

An iV-body system is described with the microscopic Hamiltonian 

N N 

H = T + V = Y,T l + £ V l3 (2.25) 

i=l i>j=l 

where for simplicity we only display two-body forces. This Hamiltonian is invariant 
under rotations, translations and reflections. We do not display nor discuss the removal 
of the centre of mass. 

At positive excitation energies, several channels may be open. In each channel, the 
particles are divided into various groups. Each such division is known as a partition. 
A given partition is denoted as a. For simplicity, we only consider here channels where 
the particles form only two subsystems: = + Nj^\ Both subsystems of partition 
a are described with an internal Hamiltonian (i = 1,2) which has the form (12.251) 
with N replaced by Nf£\ 
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A channel c is denned by specifying in addition the energy of each subsystem, i.e. 
a certain eigenvalue for each internal Hamiltonian. The exact or approximate energies 
E® and wave functions (jyg> are related by 

E^ = (^\H^). (2.26) 

Here and in what follows, subscript c has a variable symbolic meaning depending on 
the considered quantity. In Hf> = ffW, it represents the set of internal coordinates of 
subsystem i. In Eg\ it represents a set of quantum numbers. In 0W, it means both. 
When 0w i s not an exact eigenfunction of ( I2.26P remains valid from a variational 
perspective. 

To each partition a of the system into two subsystems may correspond several 
channels differing by their internal states. Each channel c has a threshold energy 

E c = + (2.27) 

defined with respect to some common reference energy. A channel is open or closed 
according to whether E c is smaller or larger than the total energy E of the system. In 



each open channel, one can define a wave number k c = J2/i c (E — E c )/h and a relative 
velocity v c = hk c / /i c , /i c being the reduced mass of partition a. In each closed channel, 
one can define a wave number k c = ^jl\i c {E c — E)/h. 

For each partition a, the relative coordinate r c = r a is the difference between 
the centre-of-mass coordinates of the subsystems. The eigenstates of each subsystem 
i = l,2 are characterized by their energy E^> and by their good quantum numbers, i.e. 
the total angular momentum li (usually called spin) and its projection M,. Under the 
time reversal operator K [50], they are assumed to transform according to the convention 

K\JM) = (-l) J ~ M \J - M). (2.28) 

A channel state is an eigenstate of the total angular momentum of the full system 
resulting from the coupling of both internal states 0c/| Ml and 4>^ 2 m 2 with a spherical 
harmonics depending on the angles Q c defining the orientation of the relative coordinate 
r c . More precisely, a channel state is represented as 



\ c ) = i'c s (2) ]/c ® y fc (n e )j"" , (2.29) 

where I c is the channel spin resulting from the coupling of I± and I2, l c is the orbital 
momentum of the relative motion in channel c, J is the total angular momentum 
quantum number of the many-body system, M is its projection, and ir is the total 
parity. Thanks to the phase factor i lc , the channel states transform under time reversal 
according to (12.281) . For simplicity, we do not display the parity quantum numbers 7r c i 
and 7r C 2 of the subsystems. They are related to the total parity by 

7T = 7T cl 7r c2 (-l)^. (2.30) 

Channel states are assumed orthogonal to each other and normed, 

(c\c') = 5 cd . (2.31) 



JM% 
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The orthogonality is not obvious for different partitions a and a'. In that case, it is 
only true asymptotically. 

Since Hamiltonian H is invariant under rotation and reflection, J, M and n are 
good quantum numbers. A partial wave of the total wave function of the system at 
energy E can be written as 

= E^-V (co) (r c ), (2.32) 

c 

where indices cq and c may represent either all quantum numbers appearing (or 
understood) in the right-hand side of (12.291) or a relevant subset of them. Projector 
A performs any antisymmetrization due to the indistinguishability of some identical 
particles. These can be electrons in atomic physics or nucleons in nuclear physics 
within the isospin formalism. The subscript (cq) recalls the entrance channel as 
explained below. This wave function is an approximate eigenstate of the full many- 
body Hamiltonian 

Mgf = J5¥gf. (2.33) 

To complete the definition of ^fffi, one must specify its asymptotic behaviour. 

We shall describe the asymptotic behaviour in terms of the collision or scattering 
matrix which is usually denoted as U in the i?-matrix context. This matrix is interesting 
physically because it is directly related to cross sections but it leads to complex radial 
wave functions. Calculations only involving real radial wave functions are also possible 
(for real potentials). They are dominantly used in atomic physics [5T] but are also 
encountered in nuclear physics [9]. The relation between both approaches is summarized 
in appendix lAl 

The asymptotic behaviour in open channels generalizing (j2.19j) is given by 

""fooW r^c C c v: 1/2 [5cc I c (k c r c ) - U CCQ O c {k c r c )}. (2.34) 

In f!2.34p . I c and O c are defined as in (12.81) and C CQ is arbitrary (see appendix |A}. The 
asymptotic form is chosen in such a way that incoming flux only occurs in the entrance 
channel Cq. Taking all possible entrance channels Cq into account, the coefficients U cco of 
the outgoing waves form the collision matrix U. It is defined for each angular momentum 
J and parity ir. Its dimension is given by the number of open channels at energy E. 
Closed channels may also appear in expansion ( I2.32p but the asymptotic behaviour of 
the corresponding radial functions is exponentially decreasing according to (I2.10p . 

For real potentials, thanks to the introduction of coefficients f~ 1//2 , current 
conservation imposes that the collision matrix is unitary, 

UU ] = U ] U = 1. (2.35) 

Because of time-reversal invariance, it is also symmetric, 

U = U T (2.36) 

where T means transposition. This property imposes the phase i lc in the definition (12.271) 
of the channel states. As shown by Huby [52], this factor is missing in some important 
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references [31 [53] and some phases must be corrected accordingly. This correction is 
not necessary if the symmetry property f !2.36|) is never used. A matrix with properties 
(I2.35P and (12.36}) can be diagonalized with an orthogonal (real) matrix S, 

SUS T = e 2 * 6 , (2.37) 

where 8 is a diagonal matrix whose elements are the eigenphases 8 n . 

Introducing (12. 32p in the Schrodinger equation (I2.33P and projecting over a channel 
wave function |c) leads to the coupled equations 

roo 

[T c + V c {r) + E C -E] u c(co) (r) + £ / W cc , (r, r')u c , {co) (r')dr' = 0. (2.38) 

/ J 



Tc=~[^-^^-). (2-39) 



They involve the kinetic-energy operators 

(f_ _ uic + 1; 

2/i c \dr 2 r 2 
The local or direct potentials are defined by 

V e (r e ) = (c\V\c) (2.40) 

where V is the total potential appearing in the many-body Hamiltonian (12.251) and the 
integration is performed over the internal coordinates of the subsystems appearing in 
channel c. The non-local potentials 

W cc ,(r,r') = (c\6(r c - r)VA5(r c , - r')\c') - V c (r)5 cd 5(r - r') (2.41) 

occur because of antisymmetrization and/or because several partitions are taken into 
account. 

Mathematically, system (12. 38ft can be written as a function of a single coordinate 
r (which becomes r' in the integrals of the non-local terms). Wave functions however 
depend on several relative coordinates r c when several partitions are taken into account. 
Channels c differ either by the nature of the subsystems or by their level of excitation. To 
simplify the presentation, we now consider that a single partition is taken into account 
or that all coordinates r c are approximated by a single one. Interactions may still be 
non local but this does not raise major problems as long as the non-local terms are 
short-ranged. 

The colliding systems have initial orientations specified by the spin projections M\ 
and M 2 in the entrance channel now denoted as c. One is looking for a solution of the 
Schrodinger equation with the asymptotic behaviour 



+(2-)- 3 ' 2 E rfW^SW (2.42) 



c'M'M' 



Several scattering amplitudes f^^M^ a PP ear - The partial wave functions (12.321) read 
with |c) representing \ {a.I\I-i)Ild r M) , 

*f c rW = £|cVW0- (2-43) 
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A stationary scattering wave function is constructed as 

Jtt II 

x (I 1 I 2 M 1 M 2 \IM)(IlM0\JM)^(^ l7T (r) (2.44) 

with M = Mi + M 2 . From the outgoing waves of the asymptotic form of (I2.44p . one 
deduces the scattering amplitudes 

/ c ( ^ 2) (^ 

K Jtt II VV 

x {IlMQ\JM){I[I' 2 M[M' 2 \I'M'){I'l'M'M - M'\JM) 

x {5 dc 5 ri 5 vl - Ufi, v ,cii)YF- M ' '(A) (2.45) 

where M' = M{ + Mg. The elastic cross section averaged over initial orientations and 
summed over final orientations reads 

da e i 1 

~oW ~ (2/ 1 + l)(2/ 2 + l) 

x E E I/c^^m^m^ + ^^^I 2 (2.46) 

MiM 2 M[M' 2 

where fc is defined in ( 12. 13ft . Inelastic or reaction cross sections are given by 

"Tsr = (2/, + m, + 1) j£. i^'( n)|2 (2 - 47) 

for c' 7^ c. The summations can be performed analytically. A long but simple calculation 
provides [53] 

d ^f-F (2/ 1 + 1 )(2/ 2 + l) ?^ )fi(C ° 3g) - (2 ' 48 ' 
where the anisotropy coefficients B\(E) are given by 

B X (E) = ^EEE E (-l) 7 ^V^+^™>Z(ZJLJ',/A) 

x Z(l'JL'J', I'\)U J d l VtCll {E) Utt' L \ cIL (E). (2.49) 

The real coefficients Z defined in [53] are modified here for consistency with the 
symmetry property (12. 36p of the collision matrix as [52] 

Z(ULJ', IX) = (-1) J+J '[(2A + 1)(2J + 1)(2J' + 1)(2Z + 1)(2L + 1)] 1/2 
/ Z L A A f I L A 1 

x V o o ojjj' j / ]' 

They verify the symmetry relation 

Z(ULJ', IX) = Z(LJ'U, IX). (2.51) 
From Bq, one deduces the integrated inelastic or reaction cross sections 

= P (2j 1+ i)k + i) P 2J+1) ¥^ fe - j,(fi)l '- (2 ' 52) 



(2.50) 
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For the elastic cross section, the summation over the orientations can also be performed 
analytically [53] and provides three contributions: nuclear and Coulomb cross sections, 
as well as an interference term. However, in practical applications, it turns out that 
definition ( 12.460 is more direct to use. 



2.4- Radiative capture 

In nuclear physics, radiative capture is an important process because of its astrophysical 
applications [51] • In this reaction, the two colliding nuclei fuse into a nucleus with 
mass m with the emission of a photon. In stars, this process often occurs at very low 
scattering energies and requires that the reaction has a positive threshold energy Q = 
(nil + m 2 — rn)c 2 . The analog reaction in atomic physics is electron capture. However, 
much more attention is paid in that field to the reversed process, namely photoionization. 
Here we shall proceed with radiative capture in a nuclear physics context but the 
formulas can be easily adapted to photodissociation or to photoionization by using 
the detailed balance. 

Radiative capture is an electromagnetic transition from a scattering state to a 
bound state. The electromagnetic aspects of this process can be studied at first order 
of perturbation theory [5J, with an outgoing scattering state ^J A . /lM2 ) as initial state 
at positive energy E and a bound state in partial wave J/vrj as final state at negative 
energy Ef. The final energy Ef is equal to — Q + E x where E x is the excitation energy 
of the final level with respect to the ground state. The radiative capture cross section 
to this state is then given by 

64vr 4 1 A + l 

ajfnf{ ] hv (2/ 1 + l)(2/ 2 + l)^[(2A + l)!!] 2 A 

x J2 \(* JfMfnf \M° x \^ MlM2) (E))\ 2 (2.53) 

M 1 M 2 Mffj, 

where the symbols aX label the electric (EA) and magnetic (MA) multipoles and the 
corresponding multipole operators are denoted as -M^ A {ji = —A, +A) [55] . The photon 
wave number is given by 

ky = {E- E f )/hc. (2.54) 

In practice, the sum over aX can usually be restricted to the dominant electric multipole 
(El, or E2 when El is forbidden) because /c 7 is small with respect to the inverse of the 
nucleus dimension. Below the Coulomb barrier, the cross section strongly depends on 
energy. To reduce the energy dependence, it is often converted into the astrophysical 
S'-factor, defined as 

Sj m (E) =Eexp(2n7 1 )a Jfnf (E), (2.55) 

where rj is the Sommerfeld parameter. 

Let us restrict the scattering wave function to a single channel. The spins ii and 
I2 are then fixed and only / and I are needed to specify the entrance channel. By using 
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expansion (12.441) in this particular case, the cross section (12. 53ft can be written as 

a Jf7Tf (E) = J2^ f (E), (2.56) 



where the partial cross section of multipolarity a A reads, in analogy with (I2.52p . 



7T 1 

& 



°>^ E) = ^ p77TT)kTT) g (2J+1) ?l^ (i; '^ - 2 ' (2 - 57) 

In (I2.57p . tiff is dimensionless and proportional to a matrix element of the 
electromagnetic operator between the final state and the initial partial wave. From 
( jgHD and f l233|) . it is given by 

2J / + i y/ 2 / 87r(A + l)^ A+1 \ 1/2 
2J+17 V^;A(2A + l)!! 2 J 

x ^(^/||^||^ (£)), (2.58) 
wz 

where the reduced matrix element is defined by 

(y J f M f*f\M* x \V{%*) = {J\Mii\J f M f ){y J t*f\\M' TX \\$> J ( J l) ). (2.59) 

For a number of reactions involving light nuclei, capture mainly occurs at distances 
where the wave functions of the colliding nuclei overlap weakly. This situation can 
be described by a simple model where the internal structure of the colliding nuclei is 
neglected and the physics of the process is modeled by a local potential V depending 
on the distance r between the centres of mass of the nuclei. In this case, the asymptotic 
form of the initial state is described by f)2.42[) where the internal states </>c/j Ml an d 4>ci 2 M 2 
reduce to spinors \I\Mi) and \I 2 M 2 ). 

The electric operators -M^ A are given in this simple model by 

Mf = eZifVYm (2.60) 
where Z^ x ^ is the effective charge 

r(EA) _ ry ( ™2 \ A . ry ( ^ > 



U? x (E,Jn^J f n f ) 



The radiative- capture cross section to a final bound state with angular momentum Jf 
can be calculated in this model. The initial scattering state with quantum numbers JM 
is defined by (12 .44ft where (12.431) is replaced by 

ty JM (r) = J2^\(hl2)njM)r^ 1 uj l (r) (2.62) 
n 

and radial functions are normalized according to (I2.34p . The normed final bound state 
with quantum numbers JfMf is assumed to be approximated by expression (12.621) with 
the normalization 



I,! 



£/ u/' (r)*dr = l. (2.63) 
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The reduced matrix element reads 

(&™\\M EX \\^{J2)) = eZ^[An(2J f + l)]" 1 ^ 

roc j 

x ^(-l) If - J Z(l f JfkJ,JA)j o u J / flf {r)r x u J IMIl) {r)dr. (2.64) 
iflfk 

In practice, (I2.64p must often be corrected empirically by multiplicative factors called 
spectroscopic factors to take an approximate account of the internal structure of the 
nuclei [56] . 



The R-matrix Theory 



17 



3. The calculable R matrix 

3.1. Introduction 

The two variants of the R matrix method mainly differ by their types of applications. 
In the calculable R matrix, the aim is to accurately solve a given Schrodinger equation 
mostly in the continuum, i.e. for positive energies. In the phenomenological R matrix, 
the goal is to parametrize scattering data; it is thus essential to know the analytical 
form of the R matrix. Of course, both variants have much in common and it is a matter 
of taste to start with one or the other. Historically, in nuclear physics, the emphasis 
has first been put on the phenomenological variant. Conversely, most applications in 
atomic physics are related to the calculable R matrix. Here we will start with a general 
formalism leading to the calculable version and then deduce the properties allowing the 
phenomenological use. 

First, we restrict ourselves to a single channel for an arbitrary partial wave. This 
assumption is often valid, and allows simple notations. We thus attempt to find 
approximate solutions of the Schrodinger equation for the relative motion of two particles 
with reduced mass /i interacting via a central potential V. At large relative distances 
r, the interaction reduces to the Coulomb interaction Vc- 

After separation of the angular part, the radial Schrodinger equation (j2.16j) for 
partial wave / can be written as 



where TJ is given by (12.391) . We are interested in bounded solutions w;(r) of (13. 1 p verifying 
condition (I2.17P at the origin. Bound-state solutions at negative energies are square 
integrable over (0, oo) and can be normed. Scattering solutions at positive energies are 
assumed to be normalized according to (j2.19j) with the scattering matrix Ui defined 
in (12 . 201) . We will essentially deal with real potentials for which the phase shifts are 
real and the scattering matrix is unitary. The generalization to complex potentials is 
straightforward. 

3.2. Definition and calculation of R matrix 

In the i?-matrix method, the configuration space is divided at the channel radius a into 
an internal region and an external region. The channel radius is chosen large enough 
so that V can be approximated by Vc in the external region at the required accuracy. 
This means that the channel radius can in principle always be increased although often 
at a cost of computational time. At each energy E, the wave function is defined by 
different expressions in these regions. In the external region, the wave function ui(r) is 



(Hi - E) Ul = 0. 
In this expression, the radial Hamiltonian Hi is defined as 



(3.1) 




(3.2) 
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approximated by the exact asymptotic expression (I2.19jj j. 

< xt (r) = Qikikr) - UiOiikr)}. (3.3) 

In the internal region, the wave function uj nt (r) is expanded over some finite basis 
involving N linearly independent functions (fj as 



iV 



wj nt (r) 



Ew(r). (3-4) 

i=i 

The functions cpj must vanish at the origin but are not necessarily orthogonal. In 
contrast with some traditional presentations of the i?-matrix theory [8], we do not 
assume that they satisfy specific boundary conditions at r = a. Various choices are 
possible, as exemplified in section HI The internal and external pieces of the radial 
functions will be connected at the boundary a by the continuity of the wave function u\ 
and of its first derivative. 

The R matrix at energy E is defined through^ 

u,(a) = Ri{E)[aui{a) - Bui(a)\. (3.5) 

A dimensionless boundary parameter B is included for later convenience. Its choice 
will be discussed later. The inverse of the R matrix is thus the difference between the 
logarithmic derivative of the radial wave function at the boundary between both regions, 
and the boundary parameter B. This matrix has dimension 1 in a single- channel case 
and is just a function of energy. In multichannel problems, the dimension of the R matrix 
is equal to the number of channels (see section [3. 101) . The principle of the method relies 
on the facts that the R matrix can be calculated from properties of the Hamiltonian in 
the internal region and that its knowledge allows determining the scattering matrix in 
the external region. 

The operator Hi is not Hermitian over the internal region (0, a). This property is 
not convenient for practical resolutions of the Schrodinger equation. This problem is 
elegantly solved with the help of the surface operator introduced by Bloch [11] 

= J «r- ) (3.6) 

The operator Hi + C(B) is Hermitian over (0, a) when B is real [TT]. Moreover it has a 
fully discrete spectrum as the self-adjoint problem is defined over a finite interval. 

The Schrodinger equation in the internal region is approximated by the 
inhomogeneous Bloch-Schrodinger equation 

(Hi + C(B) - E)u\ nt = C(B)uf\ (3.7) 

where the external solution is used in the right-hand member. The mathematical 
problem is complemented with the continuity condition 

u ) nt (a) = u1 xt (a). (3.8) 

| From now on, the Sommerfcld parameter rj is implied. 

§ As denned here, the R matrix is dimensionless. In some works, it has the dimension of a length and 
differs from the present definition by a factor a. The definition (|3.21|) of the reduced width amplitudes 
must then be modified accordingly. 
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Until now, the approximation only consists in using in the right-hand side of ( 13.71) the 
asymptotic form (j2.19|) which is known except for the value of the scattering matrix U\. 
The main advantage of the i?-matrix method is that an expansion in square-integrable 
functions can now be used in the internal region. 

Because of the Dirac function in the Bloch operator, (13. 7p and ( 13.81) are equivalent 
to the Schrodinger equation (13. ip restricted to the interval (0, a) supplemented by the 
continuity condition at r = a |llj . 

uj nt '(a) = uf^'O) (3.9) 

for any B. Hence, beyond making Hi + C(B) Hermitian, the Bloch operator enforces 
the continuity of the derivative of the wave function. The importance of this aspect 
of the Bloch operator has often been underestimated in the literature. Condition ( 13.91 ) 
needs not be imposed to the basis functions ipj since the Bloch operator will impose it 
to the physical solution u\. For historical reasons, a lot of confusion about the R matrix 
arose from the misunderstanding of this property as we shall see in section 13 .51 

Formally, the inhomogeneous Bloch-Schrodinger equation (j3.7j) can be solved with 
the Green function defined by 

{Hi + C(B) - E)G l {r, r') = 5(r - r') (3.10) 

and Gi(0,r) = 0. The solution reads 

uf lt (r)= P Gi(r,r')C(B)u1 xt (r')dr'. (3.11) 



With ( 13.61) and (I3.5|) . the R matrix is thus given by 

h 2 

Ri(E) = —Gi(a,a). (3.12) 
2fia 

The calculable i?-matrix method consists in solving the Bloch-Schrodinger equation with 
an approximate Green function expanded over a finite basis. 

To obtain a practical expression for (I3.12p . expansion (13. 4p is introduced in ( 13.71) 
and the resulting equation is projected on <fii(r), giving 

]T Cij(E, B)c 3 = —<pi(a) (au^'ia) - Buf«(a)) . (3.13) 

j—l ZfiCL 

The elements of the symmetric matrix C are defined as 

C l3 (E,B) = ( Vl \T l + C{B) + V -E\ V] ). (3.14) 

Dirac brackets correspond here to one-dimensional integrals over the variable r from 
to a. Because of the Bloch operator, the right-hand side of (13 . 1 3[) only involves values 
at r = a. 

Coefficients Cj are obtained by solving system (I3.13p . Introducing them in (13 .4p at 
r = a and comparing with (13. 5p provides the calculable R matrix 

h 2 N 

Ri(E,B) = — £ Vi (a)(C- l ) im (a). (3.15) 
-/"'' i,j=i 
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This expression is nothing but a finite-basis approximation of (13.121) . 
The wave function in the internal region is given by 



2 TV TV 



" |nt(r) = 2 f mi(E 1 B) UfXt{a) ?/ i(r) YXC-^mia). (3.16) 
We shall see in section that it does not depend on B. 

3. 3. Properties of the R matrix 

Temporarily, in this section, the basis functions (fi(r) are assumed to be orthonormal. 
The matrix of overlaps (filfj) is thus the unit matrix. Let us consider the eigenvalues 
E n i and the corresponding normalized eigenvectors v n i of matrix C(0,B), 

C(0,B)v nl = E nl v nl (3.17) 
with the orthonormality property 

VrfV n n = <W- (3.18) 
With the spectral decomposition 

TV T 

[C^BT^j:^ 1 ^-, (3.19) 

n=l h nl - h 

the R function ( 13.1 5p becomes 

TV 2 

^. B ) = Ef ! S ( 3 - 2 °) 

n=l 

with 

-"" = Wa) Ma) (3 - 21) 

and 

TV 

Mr) = ^2v n i,iLpi(r), (3.22) 
i=i 

where v n ^i is the ith component of v n i- The 7„; are known as the reduced width 
amplitudes and their squares 7^ as the reduced widths [5]. Their interpretation 
is simple. They are proportional to the value at the channel radius of variational 
approximations <p n i of the eigenfunctions of the Hermitian operator Hi + C(B). Those 
corresponding to the lowest energies thus represent approximate eigenfunctions of the 
physical problem confined over the interval (0, a) with logarithmic derivative B at r = a. 

Expression (I3.20p looks familiar to practitioners of the .R-matrix theory. It is 
however obtained here with a finite basis. The traditional expression for the R matrix 
is obtained when N tends towards infinity in a complete basis as 

= £ p-*- (3.23) 

71=1 

The energies E ni are now the exact eigenvalues of the operator Hi+C(B) and the reduced 
width amplitudes 7^ are related to the values at r = a of its exact eigenfunctions. 
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The R matrix is a real function when V and B are real. It has an infinity of real 
simple poles, bounded from below. Its derivative is always positive at regular points. 
It is a meromorphic function of the energy when the energy is considered as a complex 
variable [5]. All residues are negative and given by minus the reduced widths 7^. 



3.4- Scattering matrix and phase shifts 

Since the R matrix is known, the external function (13.31) can be introduced in relation 
(13. 8 p to determine the scattering matrix for the /th partial wave as 

l-(L l -B)R l (E,B)- [6 - M) 

In this expression, 

, OUka) 

L ' = ka 4d (3 ' 25) 

is the dimensionless logarithmic derivative of 0\ at the channel radius, L\ is the 
conjugate of L h and 

<pi = arg Ii(ka) = — arctan[Fi(ka)/Gi(ka)] (3.26) 

is the hard-sphere phase shift. Note that the same notation <fri in j5] represents the 
opposite of the hard-sphere phase shift. 

Expression (I3.24p has the striking property that it does not depend on the boundary 
parameter B, independently of the size of the basis. Indeed, with the matrix relation 
(IB4I) in Appendix [B], one deduces from (13. 151) and (13. 14f) 

= Rdks) + B (3 - 27) 

for any B, real or complex. Expression (13.271) means that the logarithmic derivative of 
the internal wave function at the boundary is independent of B. Introducing relation 
(I3.27P into (I3.24p shows that any B value leads to the same scattering matrix as for 
-6 = 0. Equation (I3.27P is well known in i?-matrix theory (see equation (IV. 2. 5a) of [5]). 
It is also valid for the phenomenological R matrix with a finite number of poles [H]. 
However its validity for the approximation (I3.24p for any basis size [2H1 EZj is sometimes 
overlooked. 

Like the scattering matrix Ui and the external wave function ttf xt (r), the internal 
function (I3.16P does not depend on the choice for B. Indeed, with the help of relation 
flB3p . one easily shows that, for any B, it is equal to the similar expression where B is 
replaced by zero. 

For a better physical interpretation of the results which is important in applications, 
it is convenient to introduce some definitions. To this end, Li is separated into its real 
and imaginary parts as 

Li = Si + iP L . (3.28) 
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The real part Si and imaginary part Pi of Li are called the shift and penetration factors, 
respectively. They depend on energy and on the channel radius. The penetration factor 
can be written with the Wronskian relation liO\ — I[0\ = 2i as 

/ 7— 1\ kct Jvd , , 

' ( )= \Oi{ka)\* = Fi{kaY + Gi{kay (3 - 29) 
It is always positive and increasing [5]. The shift factor reads 

Si{E) = Pi{E)[Fi{ka)Fl{ka) + G^G'^ka)]. (3.30) 

It is always negative for rj > [5j. Although we do not know a proof that Si is always 
increasing in the same case, we could not find numerically a counterexample. As shown 
below, none of these properties is valid in the attractive case. 
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Figure 1. Penetration factors Pi{E) (upper panel) and shift factors Si(E) (lower 
panel) in the neutral case (rj = 0) as a function of E in units of fr 2 /2/ia 2 . 



In the neutral case (77 = 0), the penetration factors have simple analytical 
expressions such as 

P (E) = ka, Pi(£) = 1 j^ )2 , ... (3.31) 

They do not vary very fast with energy (see Fig. [1]). This figure is universal, i.e., 
independent of the collision. Notice that the derivative of Pq with respect to energy is 
infinite at the origin. This property leads to the special behaviour of neutron scattering 
in the s wave. Penetration factors decrease with the orbital momentum / as expected 
from the occurrence of an increasing centrifugal barrier. The shift factors read 

s (£) = o, ^) = -rrkr ■■• (3 - 32) 
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They vary smoothly with energy, starting from the integer values —I (see Fig. [I]). This 
weak energy dependence is the origin of the Thomas approximation j5] where the shift 
factor is assumed to vary linearly in a limited energy range. 
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Figure 2. Penetration factors Pi(E) (upper panel) and shift factors Si(E) (lower 
panel) in the repulsive charged case for a = as a function of E in units of h 2 /2fia 2 . 



The penetration factors are very different according to whether both particles are 
charged or not. In the repulsive charged case, the energy dependence of the penetration 
factors is much stronger (notice the logarithmic scale in Fig. [2]). Here we have to choose 
the strength of the Coulomb interaction. Figure |2] corresponds to a channel radius a 
equal to the Bohr radius a#. The strong dependence at low energies is due to the 
difficulty of penetrating a Coulomb barrier when the scattering energy becomes much 
smaller than the Coulomb barrier. Beyond 1 = 1, increasing from I to I + 1 decreases 
the penetration factors by more than an order of magnitude. In contrast, shift factors 
are much more similar in all cases. In the repulsive charged case, the shift factors are 
quite similar to those of the neutral case, except for I = (see Fig. 

In the attractive charged case, the energy dependence of the penetration factors is 
displayed in Fig. [3] for a = as- It is rather similar to the neutral case except at low 
energies where it starts from a finite value at energy zero. Shift factors displayed in 
Fig. [3] resemble other cases. Notice however that Sq is positive, and decreasing. 

With definition (I3.28p . the collision matrix (13.241) becomes 

i5i _ , ; i - {Si - m + 

i-{s l -B)R l -im' [ } 
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Figure 3. Penetration factors Pi(E) (upper panel) and shift factors Si(E) (lower 
panel) in the attractive charged case for a = as as a function of E in units of h 2 /2^ia 2 . 



u) nt (r) 



One obtains an explicit expression for the phase shift, 

P R 

5i = (j)i + arctan — - 1 . (3.34) 

1 — (bi — D)Ki 

The internal wave function (I3.16P can be rewritten as 

h 2 1 (Jen PA 1 / 2 N N 

TTa 6 " ^rOiffi E ^(0 UC~%m(a), (3.35) 
fJ-a \L — {Li — n)Ki\ j=1 i=1 

which explicitly shows the phase and modulus of wj nt (up to a global sign). For an 
orthonormal basis, this expression can be rewritten using (I3.19P and (13. 22ft as 

<{r) = £>-frq , {k ( f^ m] £ ^fl. (3.36) 

3. 5. On the basis and boundary parameter choices 

Considerable confusion exists in the literature about the properties that basis states 
(Pi should have. Improper choices have led to the introduction of corrections and to 
attempts to use the boundary parameter B to correct drawbacks of the basis. However 
we have just seen that the results are independent of the choice of B. It is thus 
worthwhile to devote this section to a clarification of this issue that has sometimes led 
to an undeserved reputation of poor convergence for the calculable i?-matrix method. 

In their seminal paper, Wigner and Eisenbud wanted to provide a phenomenological 
description of resonances [3]. They did not intend to propose a technique of resolution. 
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To reach their goal they assume that the basis functions all satisfy (for B = 0) the 
boundary conditions <Pj(0) = and 

aip'^a) - B(pj(a) = 0. (3.37) 

This procedure leads to R matrix (I3.23p . When used as a technique of resolution, the 
finite-basis R matrix (I3.15P or ( I3.2(jp obtained with this procedure does not converge 
uniformly. The reason is simple. The first derivative of the wave function C I3 .161) suffers 
from a discontinuity at r = a [58] . The limit of wj nt when r tends towards a to the left 
is equal to u^ xt '(a) but not to Mj nt '(a), 

lim u) nt '(r) = ^ xt '(a) ^ ^'(a). (3.38) 

r— >a~ 

For example, if B = 0, (f'j(a) vanishes for all j values and one readily sees that u\ nt '(a) = 
at all energies. This property has unfavourable consequences on the convergence of 
numerical methods when the basis is truncated since the logarithmic derivative of the 
external solution depends on the phase shift (and thus on energy) and can not be 
matched with the internal solution (see Figures |5] and [8] in section H]). Buttle [27] has 
proposed a correction to the i?-matrix truncation. His idea is to replace the truncated 
part by an analytical approximation, i.e., in practice, by the same expression for the 
zero potential. Although this correction improves the phase shifts, it does not really 
solve the problem because it does not improve the wave functions. 

This problem received a solution with the works of Lane and Robson [281 [291 159] . 
Their method was successfully applied in nuclear physics where traditional basis 
functions do not satisfy (I3.37P and, on the contrary, display a variety of behaviours at 
the channel radius. With oscillator basis functions, accurate results for neutron-nucleus 
scattering could be obtained [5H ED]- At the same time, a microscopic extension of the 
R matrix using a fully antisymmetrized two-centre harmonic-oscillator model provided 
accurate phase shifts for collisions between light nuclei with few basis states EE] (see 
section 16721) . The success of these calculations relies on the fact that the Bloch operator 
makes condition f!3.37j) unnecessary. Since the results do not depend on B, the choice 
B — was used. A general though economical method for solving coupled-channel 
problems is described in [H]. 

The negative role of condition ( I3.37P remained long unnoticed in atomic physics 
where in many cases basis states were imposed to satisfy (13.371) [6T1 El EH E2J [63] . 
In the literature, the basis functions are often assumed to be solutions of a Sturm- 
Liouville problem satisfying (I3.37p . As reviewed in [9], various solutions to the purported 
convergence problems of the calculable R matrix have been proposed, such as using two 
different sizes for the internal region. Another type of solution proposed in [M] requires 
a basis with a boundary condition depending on the eigenvalue E n \. In 1983, Greene 
[65. 66 realized in the context of atomic physics that, in place of the traditional choice 
(I3.37P of a common boundary condition to all basis functions at r = a, a variety of 
values for the logarithmic derivatives of the basis functions should be far more efficient 
for accurate calculations. He also proposed to optimize the boundary parameter B 
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to have a better connection between the internal and external logarithmic derivatives. 
With his variational principle (66] applied to potential scattering, B is given at each 
energy E by the generalized eigenvalue problem 

C(E,B)c = (3.39) 

in the present notation. One readily sees [H7] that the only eigenvalue B is 1/R(E) 
which, as the logarithmic derivative of the internal wave function, is obviously optimal. 
With this choice, the basis functions allow a perfect matching but at the cost of a new 
calculation of the basis at each energy. As we have shown with (13.271) . this complication 
is unnecessary since all physical results are independent of B. 

The first accurate calculations in atomic physics with a basis that fully ignores 
condition (I3.37H can be found in [6E1 ESS, EHJ [70] . New developments making use of B 
splines now avoid this condition (TTJ [721 EH]- So let us emphasize that the calculable 
R matrix does converge accurately when the basis functions are a well chosen part of 
a complete set displaying a variety of logarithmic derivatives at the channel radius a. 
The reason of its accuracy is that the Bloch operator imposes a good matching at the 
boundary [ID]. The choice of a boundary parameter B is irrelevant and the Buttle 
correction is not necessary because the good matching of the internal and external wave 
functions allows an accurate determination of the phase shift. Practice has shown that 
this simple solution allows a much smoother connection with the external wave function 
[4"U| |4"T] . In opposition to the traditional presentation, there is thus no need for a special 
assumption about the behaviour of the basis functions at the boundary. 



3.6. Resonances 

Resonances can be studied in various ways. Each of them may be useful, either in 
calculable applications (section B~8j) or in phenomenological applications (section [5]). 
In a first approach, the boundary parameter can be chosen as [TH [2H] 

B = U. (3.40) 

It thus depends on energy. This complex value leads to a complex function Ri(E,L{) 
which is not an R matrix in the strict sense since B is not a real constant. It is the 
function introduced by Kapur and Peierls [I]. Nevertheless it is also given by expression 
f)3.15p . Equation (13.241) then takes the simpler form [2S1 EH] 

Ui = e 2 ^ [1 + (Li - Li)IU(E, U)]. (3.41) 

This expression is also valid for complex k values if <pi is defined as the phase of 
Ii(ka). Since (13.41 j) has no denominator, a direct relation appears between a pole of 
the scattering matrix, i.e. a resonance energy, and a pole of the complex R matrix. 
This relation is however valid only when Li is calculated at the resonance energy. This 
means that the scattering and R matrices have only one common pole and only at 
specific energies. Determining S-matrix poles in this way thus requires some iterative 
procedure. The choice (I3.40p for the boundary parameter can be used to analyze the 



The R-matrix Theory 



27 



mathematical nature of resonances in the complex plane in coupled-channel cases [71] . 
However, the same results can also be interpreted with the traditional R matrix involving 
only real parameters [75]. Let us return to real energies for the other approaches. We 
choose B = to simplify the presentation. 

Another definition of a resonance energy Er is that it corresponds to the value 7r/2 
of the resonant part 5i — 4>i of the phase shift. From ( I3.34p . it is therefore defined by 
the equation 

1 - S^Eh^Er) = 0. (3.42) 

In general this equation must be solved numerically. To define the resonance width, let 
us consider the collision matrix f!3.33j) for energies close to Er. A Taylor expansion of 
Si(E)Ri(E) for E ~ Er provides the Breit-Wigner approximation 

In this expression the (energy-dependent) width of the resonance is given by 

T{E) ~ [d(S lRl )/dE] E=ER (3 - 44) 
Because of the shift of E R with respect to a pole (see below), Ri(E) can be supposed to 
vary slowly near a narrow resonance (see section [5]). The total width then reads 

T(E) = 2 7 2 P,(£) = 3Wr(E R ), (3.45) 

where T(Er) is the width at the resonance energy. The reduced width 7 2 defined by 
(I3.45P is given by 

7 2 = R l {E R )/[d{S l R l )/dE} E=ER . (3.46) 

Let us mention that (13.421) may have solutions which do not correspond to physical 
states. The width of a physical resonance should be small enough to make its lifetime 
longer that the typical collision time. 

In another way of studying a resonance, let us consider an energy very close to a 
pole E n i of the R matrix. If all terms with n' ^ n can be neglected, the R matrix is 
approximated as Ri(E, 0) ps ^J{E n i — E). A simple calculation provides 

This expression resembles the Breit-Wigner form of the phase shift 

5f w ^ fa + arctan jT^l • (3-48) 

Er - E 

By comparison, one defines the resonance energy 

E R = E rd -^S l {E R ) (3.49) 
and the formal width 

T(E)=2 1 2 nl P l (E). (3.50) 
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The resonance energy is defined by an implicit equation which can be solved 
approximately (see section 15. 2p . The width is an energy-dependent quantity whose 
asymmetric shape depends on the behaviour of Pi. Its relation with a measured width 
is also discussed in section 15.21 

With (ET2B . ( jS3SD and (ET50j) . the internal wave function (I375S|> can be 
approximated at the vicinity of a resonance by 

1 1/2 

<M0- (3-51) 



It is thus proportional to an approximate eigenfunction ( I3.22p of H L + £(5) with a 
proportionality factor exhibiting the usual Lorentzian energy dependence of a resonance. 
Equation (13.5 1 p is at the basis of the so-called bound-state approximations where the 
resonance is described by a square-integrable wave function [76] . 

3. 7. Bound states 

The .R-matrix formalism can be extended to bound states (Eb < 0) [37] . In that case 
the external wave function uf^ is given by 

uT\r) = QWi(2K B r), (3.52) 

where Wi(x) is a shorthand notation for the Whittaker function (12.101) and where Kb 
and rjB are the wave number and Sommerfeld parameter, respectively, of the bound 
state. In (I3.52p . C\ is the asymptotic normalization constant (ANC) which determines 
the amplitude of the wave function at large distances. This quantity plays an important 
role in some nuclear reactions of astrophysical interest [TTJ (see section loTTj) . 

To determine Eb, & convenient choice for the boundary parameter B in ( 13.61) is 

B = U(E B ) = St(E B ) = 2« B a ^^4 , (3.53) 

Wi{2K B a) 

because it suppresses the right-hand side of the Bloch-Schrodinger equation (13. Tj) . Since 
the wave function is real, Li is real and identical to the shift factor. The penetration 
factor Pi vanishes. The internal wave function is expanded over a basis, as in (13. 4p . 
Projecting the Bloch-Schrodinger equation (13 ,7p over a basis function ipi provides for 
z = l,JV, 

N 

£ (<pi\Ti + C(Li) + V - E B \iPi) Cj = 0. (3.54) 

3=1 

This system of equations is similar to a standard eigenvalue problem, but parameter Li 
depends on energy Eb- In practice one starts from Li = and iterates until energy Eb 
has converged. At convergence, the Cj can be calculated by solving the system. If uf 11 
is normed according to (I3.18P , the square of the norm of the wave function is given by 

roc 

Ni = l + (Ci) 2 / (Wi(2K B r)) 2 dr. (3.55) 
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This expression can be rewritten [17] as 



Ni = 1 + 7. 



ni 



dSi 
dE 



(3.56) 



E=E f 



where 7„; is the reduced width amplitude of the bound state. 

The internal wave function is given by (13.161) multiplied by N l 



-1/2 



u) nt (r) 



N, 



-1/2 



N 



(3.57) 



3=1 

where coefficients v n z,j correspond for an orthogonal basis to the eigenvector v n i of 
C(0, L;) at energy E n i = Eb in ( 13 . 1 7f) . Using (13. 22ft . it can be rewritten as 



u) nt (r) 



From 



and (13321 . the ANC is given by 



or with (jgjnP by 

Q = (2 l 2a/h 2 N l ) 1 / 2 lnl /W l (2K B a) 



(3.58) 
(3.59) 
(3.60) 



It should be independent of radius a. This relation between the ANC and a reduced 
width amplitude which corresponds to a vanishing width can be useful for the 
phenomenological R matrix. A similar formalism can be applied to resonances [79] . 
and provides widths as well as energies of resonances. 



3.8. Capture cross sections 

The determination of capture cross sections requires the calculation of matrix elements 
of the electromagnetic multipole operators M. aX between an initial scattering state and 
a final bound state. One can take into account the division of the configuration space 
in the general case [37] but the principle of the calculation can be explained more easily 
in the simple potential model using equations (I2.57p . (I2.58P and (I2.64p . 

According to the i?-matrix framework, the radial matrix element between an initial 
scattering wave function Ui(r) = uf fl ^j^(r) and a final bound-state wave function 

Uf(r) = u J / fls (r) appearing in (12. 64p can be written for an electric multipole as 

/•oo pa poq 

/ u f r x Ui dr = / ufr x ufdr + / uf'r^dr. (3.61) 

JO JO J a 

The internal matrix element is given by 

f ufr^dr = J2 cwa, k P Vk>r x <p k dr, (3.62) 

where coefficients cy/j and are related to the final and initial wave functions, 
respectively. Notice that the bases could be different for both states. In the external 
region, we have 

/ uf t r x ^ t dr = QC f W lf (2K B r)r x (I k (kr) - U^ n O k (kr))dr. (3.63) 

J a J a 
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As for elastic scattering, the total matrix element ( 13. 6 1 1) should not depend on the 
channel radius a, whereas each contribution does depend on a. For transitions to weakly 
bound states, the external term is dominant [76] since the Whittaker function slowly 
decreases at large distances. Neglecting completely the internal contribution leads to 
the "external-capture" model [5U]. A typical example is the 7 Be(p,7) 8 B reaction where 
the ground state is bound by 137 keV only [77] . On the contrary, resonant reactions, or 
transitions to deeply bound states provide a dominant internal term. 

3.9. Propagation methods 

For long-range potentials, the channel radius may become very large. This may induce a 
prohibitive size for the .R-matrix basis and lead to long computation times. This problem 
was first met in atomic physics because of the long tail of polarization potentials and 
led to the development of propagation methods [321 ED E2]- This situation can also 
occur in nuclear physics, for example, in three-body systems [83] or in coupled-channel 
calculations |84j . 

Different methods have been proposed to address this problem. The basic idea is, 
either to propagate the wave function or the .R-matrix over several intervals on which 
the basis size remains reasonable, or to determine "distorted" Coulomb functions valid 
at distances shorter than the channel radius. We present here a propagation method 
directly derived for the .R-matrix (see for example [82]). 

The idea of propagation methods [32] is to divide the internal region (0, a) in N s 
subregions (a a -i,a a ) for a = 1,N S , with oq = and ajv s = a. The intermediate radii 
a a can be chosen equidistant, but this is not mandatory. The width of the intervals 
and the basis size in each interval can also depend on the number of oscillations of the 
wave function [HH E2]- With small intervals, the .R-matrix bases remain limited, but 
the number N s of repetitions of the calculation may be large. Approximations of the 
potential can often be employed in some intervals. In some cases, the size of the intervals 
is chosen small enough so that the potential may be considered as constant which allows 
an analytical propagation [32] . 

We briefly present here the principle of the propagation technique for potential 
scattering with a basis in each interval. We refer to references [HU E2] for a multichannel 
extension. A Bloch operator C a is now defined at each boundary a a as 

+ 2 j 

£ a = —5(r-a a )—, a = l,N a . (3.64) 

The Bloch-Schrodinger equation is replaced by a set of equations 

(Hi + £ Q - £ Q _i - E)ui ta = (£ Q - £ Q _i)u ZiQ , a = 1, N s (3.65) 
with £ = 0. The boundary conditions are ^i,i(0) = and 

Ul,a( a a) = Ul t a+l( a a), (X = 1, N s (3.66) 

with it^iVs+i = uf*. In each interval, the wave function is expanded over a set of N a 
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basis functions 

JVo 

^»=E C MM- (3-67) 
i=i 

Equations (13.651) can be solved with approximate Green functions as in section 13.21 By 
projecting (j3.65|) on one of the basis functions, the solutions can be approximated as 

«»,«(r) = ECCa^^lA, - C a ^\u ha )^{rl (3.68) 
ii' 

where the symmetric matrix C a is defined in each interval by 

C a>ii , = {tf\H % + C a - C a ^ - E\(p$). (3.69) 

The Dirac notation represents an integration limited to the range (a a _i, a a ). 

For each interval, (13.68)) can be used to determine a relation between values of the 
wave function and its first derivative at a Q _i and a a , 

ui, a (a a -i) = fci u i,a( a a) - 7^ii«j }a (a a _i), (3.70) 
ui, a (a a ) = TZ% u' lo ,(a a ) - TZ^u'^aa-t), (3.71) 
where various values of the approximate Green functions TZppi are defined as 

ifo = % E ^W(^VfM ( 3 - 72 ) 

with = 0, 1. 

An R matrix can be defined at each boundary with an extension of (13. 5p as 

Ui, a (a a ) = a a R(a a )u' la (a a ), a = 1,N S . (3.73) 

Equations (13.701) . (13.711) and (13. 73[) provide relationships between R matrices at a a _i 
and a a (a = 2, N s ), 

a Q _ 1 i?(a Q _ 1 ) = -n a u + KUK - a a R(a a )]- l n^ (3.74) 
a a R(a a ) = 11^ - n^[n^ + a Q _ lJ R(a Q _ 1 )]~ 1 ^ . (3.75) 

The latter equation provides an outwards propagation [R(a a ) from R(a a -i)] and the 
former provides a backwards propagation [R(a a _i) from R(a a )]. This technique is 
quite efficient in multichannel calculations. Numerically the main part of the i?-matrix 
calculation arises in the inversion of matrices C a [see f)3.72p ]. It may save computer 
time to diagonalize them when many energies are needed. The size of these matrices is 
given by the number of basis functions times the number of channels. If the calculation 
involves many channels, reducing the number of basis functions may lead to a significant 
reduction of the computer times. 

The collision matrix is obtained from R(a^ s ) with (13.24)) . The external wave 
function uf^ is thus known. By starting from the last interval, the wave function is 
determined by its coefficients in each interval. If the £ a _i terms are simplified in ( 13.65)) . 
one obtains the system 

E [CW + (<^|A*-lK')] Cj, = (^\C a \ui, a ) (3.76) 
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where f l3.T3f) can be used to eliminate ui, a '(a a ) in the right-hand side. The matrix in 
this system can be inverted with (IB2p . Starting with a — N s and going backwards one 
obtains the coefficients in all intervals. 

Let us briefly discuss other techniques dealing with long-range potentials. In [83] , 
in the framework of three-body continuum states, we propagate the wave functions from 
a a _i to a a by using the Numerov algorithm. This approach avoids the choice of a basis, 
but requires longer computer times. It is also difficult to apply to non-local potentials. 

The Light- Walker propagation consists in approximating the potential by a constant 
in small enough intervals [32]. In [M], this method is improved by considering linear 
approximations of the potential. In the method suggested by Gailitis [85] , the Coulomb 
wave functions are modified by 1/r expansions, which can be used at short distances 
(see also [M1IE2]). 



3.10. Extension to multichannel collisions 

Until now the presentation was, for the sake of clarity and simplicity, limited to single- 
channel calculations which also neglect the spins of the colliding particles. This does not 
affect the general properties of the i?-matrix theory. However many problems require a 
multichannel approach. 

In a many-body problem, the Schrodinger equation (I2.33j) involves Hamiltonian 
f)2.25p . The total wave function \fr JMn of the system, with total angular momentum J 
and parity it, is expanded over a set of channel functions (12.291) . denoted as |c). If we 
assume that a single relative coordinate appears in the problem or that we approximate 
all relative coordinates by a single one, the wave function is given by (12.431) . The 
Schrodinger equation is replaced by a set of differential equations 

J2 [(T c + E C - E)b cd + Vj\ u d = 0, (3.77) 

c' 

where, as before, T c includes the centrifugal term. This system is a particular case or a 
local approximation of (I2.38p . 

A typical example is given by coupled-channel calculations where channel functions 
are defined by (12.291) . This situation occurs, for example, in coupled-channel calculations 
with a discretized continuum (CDCC) [87], see section 16.31 Another example 
corresponds to three-body scattering, where the channel functions |c) contain various 
quantum numbers defined in the hyperspherical formalism [SB] , see section 16.41 

In all cases the problem is first to determine the potentials V cc > and then to solve 
(I3.77P for positive energies. Here we essentially deal with the second step. We thus 
assume that the potentials are known and that they present the asymptotic behaviour 

V cd — ► ^^ cc , (3.78) 

r— >oo y 

In these conditions the asymptotic form of the radial wave functions u c (r) is given by 
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( 12.341) and the radial wave functions in the external region at energy E are defined as 

cxt / \ _ \ v^ 2 (l c (k c r)6 CC0 - U CCa O c (k c r)) for E c < E 
Uc(C0)[r) \ A CC0 W_^ l+ ,(2n c r) for E c > E, ^ 

where Cq is the entrance channel (E co < E). 

The Bloch operator (13.61) is defined in the multichannel formalism as 

£ = I»A<c|, C=^r-a){±-^\ (3.80) 



where coefficients B c are chosen as zero or as in (I3.53P for open and closed channels, 
respectively. Notice that these coefficients then depend on energy for closed channels. 
This choice may be inefficient in some variants of the R matrix [59] but is convenient in 
(I3.88P and f !3.89|) below. The Bloch-Schrodinger equation is given by 

[(T c + £ C + E C - E)8 CC , + V cc ] uf{r) = £ c < xt . (3.81) 

c' 

The internal parts of the radial wave functions are expanded over a basis (fij(r), 

N 

ffl = Ew(r). (3.82) 

3=1 

The determination of the R matrix and of the collision matrix U are direct extensions 
of the formalism developed in sections 13.21 and 13.41 
The R matrix is defined as 

u c{a) = 5Z(a*c/ ' Vc>) 1,2 Rcc'[au' c ,(a) - B d u c ,(aj\. (3.83) 

d 

Matrix R is symmetric, with elements given by 

h 2 N 

Rcc'(E) = Yl <Pi{a){C-% tC/m ,(a), (3.84) 

where 

Cd >c /i' = (ifi\T c + C c + E c - E\(pi>)5 cc > + {tpi\V cc >\(pi>) . (3.85) 



Like in section 13. 3[ the spectral decomposition of the symmetric matrix C for an 
orthonormal basis provides the canonical form of the multichannel R matrix as 

RAE) = Y,¥% (3-86) 

where the real poles E n are the eigenvalues of C and the reduced-width amplitude of 
pole E n in channel c is expressed as a function of the components of the corresponding 
normed eigenvector v n as 

/ h 2 V 1/2 N 

lnc= \2iTa) Tl^^iia)- (3.87) 

The number of terms in the sum (13.861) is given by the product of the number of channels 
by N. 
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The collision matrix is obtained with 

u = z~ x z*, 

where an element of matrix Z reads 

-1/2 



Zed = {kef a) 1 O c {k c a)5 cc > - k d aR cc >0' cl {k c >a] 



(3.88) 



(3.89) 



For complex potentials as encountered in the optical model, expression (13.881) must be 
modified into 

U = Z~JZ U (3.90) 

where Zq = Z and Zj ^ Z* Q is given by a similar expression with outgoing functions 
replaced by incoming ones. Matrix U is then not unitary. 

The dimension of the R matrix is equal to the number of channels included in 
the calculation; it does not depend on energy. On the contrary the dimension of the 
collision matrix U is given by the number of open channels and may vary with energy. 
In (I3.89p . only open channels contribute thanks to the choice of boundary parameters 
B c . The R matrix can be modified by eliminating the closed channels. When not all 
channels are open, let us denote open channels by c and closed channels by c. Equation 
()3.84j) remains valid with matrix C replaced by a smaller open-channel matrix C° with 
elements [89| IB] 

C°i,cfi< = Cd,c'i' - ( l Pi\ V ce\ l Pj){C 1 )cj,&j'{ l Pj'\Vgc'\fi'), (3.91) 



where C is the restriction of the full matrix to closed channels and c, c' are the 
corresponding indices. 

The internal components of the wave function are given for both open and closed 
channels by 

N 

I' c ,{k c ,a)5 c > C0 -U c > Co O' c ,(k c/ a)] ^ V 9 i( r )(C~ 1 )ci, c vV 9 i'( a )) ( 3 - 92 ) 



<t){r 



^ h 2 k c > 

~f 2 He' a 



t,t'=i 



where the sum over d runs over open channels only. The coefficients of the external 
part of the wave function for closed channels read 



A, 



Vc,l+ 



(2« g a)] [ Y J (Hck c '/h) 1/2 aR c - c > r c ,{k d a)5 c i (X) - U c > C0 O' c ,{k c 'a) 



(3.93) 



where the sum over d also runs over open channels only. 
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4. Applications of the calculable R matrix 

4-1. Conditions of the calculations 

In this section we apply the calculable .R-matrix method to the scattering by a potential. 
Our goal is not to fit experimental data, but to illustrate the method under different 
conditions, and with different types of basis wave functions. In this approach, both 
colliding particles are assumed to be structureless, and to interact through a potential 
[see (13.2j) ]. In general, this potential V(r) involves a local term U(r) and a non-local 
term W(r,r') such that 

;>oo 

V(r)ui(r) = U(r)ui(r) + W(r,r')ui(r')dr' . (4.1) 

J 

In the following, unless specified otherwise, only the local term is included. In nuclear 
physics applications, the local potential U(r) usually contains a nuclear and a Coulomb 
contributions, denoted as V/v(r) and Vc(r), respectively. The potential may be I 
dependent, by adapting its parameters to the I value. The Schrodinger equation 
associated with potential (14.1 j) can be solved exactly (at the computer precision) with 
the Numerov method [901 19T1 192] . The exact phase shifts and wave functions will be of 
course compared with the i?-matrix calculations. 

In nuclear-physics applications, the reduced mass is expressed in terms of the 
nucleon mass as 

where A\ and A2 are the nucleon numbers. The calculations are performed with 
h 2 /2m N = 20.736 MeV.fm 2 and e 2 = 1.44 MeV.fm. 

To cover a broad variety of applications, we have first selected three typical systems: 

(i) The 12 C+p system (I = 0) which presents a narrow resonance (r = 37 keV) at 
E = 0.42 MeV. 

(ii) The a + a system (/ = 4) which presents a broad resonance (r = 3.5 MeV) near 
E = 11.3 MeV. 

(iii) The a+ 3 He system (/ = 0) which is non resonant. 

These cases are typical examples of nuclear-physics applications, and will be treated 
by different types of basis functions. Here, we only consider real potentials. This is 
consistent with the low-energy regime where the i?-matrix method is well adapted. The 
extension to complex potentials is however trivial. 

Then, the calculation of e~— H phase shifts will provide a typical example of an 
atomic-physics application with a non-local potential. Atomic units will be used. 

4-2. Basis functions 

We consider different families of basis functions <Pi(r) [i = 1, . . . , N, see (13. 4D ] vanishing 
at r = 0, commonly used in the literature. The corresponding matrix elements are given 
in Appendix C. 
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(i) Sine functions 

The basis functions are given by 

^(r)=sin— (i-1/2). (4.3) 
a 

This choice seems natural since it can simulate the oscillating behaviour of the wave 
function near the channel radius. However we will show that these basis functions 
are not suitable for accurate calculations since the derivative satisfies ( I3.37P (see 
section I3~5l) . i.e. it vanishes at r = a for each i value, 

(a) = 0. (4.4) 

This property means that any finite combination of these basis functions will present 
a zero derivative at r = a. As discussed before, the matching between the internal 
and external solutions is therefore expected to be poor, and the i?-matrix phase 
shifts rather inaccurate. 

(ii) Gaussian functions 

The basis functions have a Gaussian dependence jl3] with different size parameters 

<p i (r)=r l+1 exp(-(r/b l ) 2 ), (4.5) 

where the 6,- are chosen clS 8b geometric progression 

6 i = 6 1 4" 1 , (4.6) 

and are therefore determined by a set of 2 parameters (b\, xq). Gaussian-type basis 
functions (j4.5j) allow an analytical calculation of the overlap and kinetic-energy 
matrix elements. For Gaussian potentials, the matrix elements are analytical as 
well. 

(iii) Lagrange functions 

The Lagrange functions jl2] are defined in the (0, a) interval as 

W (r) = (-1)*+' (JLy^—, ^/*- 1) 
\aXiJ v r — aXi 

where P/v is the Legendre polynomial of order N, and Xi are the zeros of 

P N (2xi - 1) = 0. (4.8) 

In (14. 7p . the factor (r/axi) n is aimed at regularizing the basis function at the origin. 
For two-body calculations defined in the (0, a) range, we use n — 1, which ensures 
that the wave function vanishes at the origin. 
The basis functions satisfy the Lagrange conditions 

(Pi(axj) = (aAi)- 1/2 c^-, (4.9) 

where Aj is the weight of the Gauss-Legendre quadrature corresponding to the (0, 1) 
interval. This basis is exactly equivalent to the Legendre basis defined by 

< Pi {r)=rP i - 1 (2r/a-l). (4.10) 
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However, if the matrix elements with basis functions (14. 7\\ are computed at 
the Gauss approximation of order N, consistent with the N mesh points, their 
calculation is strongly simplified. At this approximation, the overlap is given by 

( { Pi\ { Pj)=l ¥i(r)Vj(r)dr « 8 i5 , (4.11) 

J 

and the local-potential matrix is diagonal with elements just given by the value of 
the potential at the mesh points, 

{ipi\U\ Vj ) = f ip^U^ip^dr « UiaxJSij. (4.12) 

J 

This calculation is easily extended to non-local potentials W(r, r') 

((fi\W\Lpj) ~ a^J Xi\jW(aXi, axj). (4-13) 

These matrix elements do not need any evaluation of integral. The kinetic energy 
term is given by analytical expressions [93] (see Appendix C). This method has 
been shown to be quite efficient and accurate in various fields for bound states 
[9"3] as well as for scattering states [121 EH]- In the present applications, the Gauss 
approximation will be used systematically, except in the e~— H scattering, where 
this approximation will be also tested by a Simpson quadrature for the matrix 
elements. 



4-3. Application to a narrow resonance: 12 C+p 

The 12 C+p system is known to present a narrow resonance (J n = l/2 + , Er = 0.42 
MeV, T = 37 keV [M]) at low energies. As mentioned before, our aim is not to fit these 
properties accurately, but to compare the exact solutions of the Schrodinger equation 
(13. lj) with the .R-matrix approach. For the nuclear and Coulomb potentials we choose 

V N (r) = -73.8exp(-(r/2.70) 2 ), 

V c (r) = 6e 2 /r. (4.14) 

In what follows, the units in the nuclear potentials are fm and MeV for lengths and 
energies, respectively. This potential reproduces the resonant behaviour of the s phase 
shifts at 0.42 MeV. To simulate the Pauli principle, a nucleus-nucleus potential may 
contain additional (unphysical) bound states [HS1 EE]- These Pauli forbidden states 
show up in microscopic calculations, where the internal structure of the colliding nuclei 
is taken into account, but their effect can be partly simulated in non-microscopic theories 
by additional bound states in the potential (see [95j |96] for details). For the 12 C+p 
system, the I = potential contains one forbidden state. 

In Fig. HI we present the exact phase shifts and .R-matrix calculations with Lagrange 
(a) and Sine (b) functions. The channel radius is chosen as a = 8 fm, where the 
nuclear interaction is negligible. Fig. H] (a) shows that with the Lagrange functions, the 
convergence is reached with N > 10. We show an example with N = 7, where N is aimed 
at being not large enough. Fig. [4] (b) illustrates the Sine functions, poorly adapted to a 
good matching since the left derivative of the wave function at r = a is zero [see ( 14.41) ] . 
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Figure 4. 12 C+p i?-matrix phase shifts (in degrees) for different bases and conditions 
(1 = 0). (a) Lagrange functions at a = 8 fm (the exact results are superimposed to the 
N = 10 curve), (b) Sine functions at a = 8 fm, (c) Convergence as a function of a for 
N = 15; the Lagrange and Gaussian results are superimposed (a = 7 fm corresponds 
to the exact results). 

Even N = 20 is far from the exact calculation. In Fig. |4] (c), we present the convergence 
as a function of the channel radius with the Lagrange functions (N = 15 is fixed). 
Results obtained with Gaussian functions (with bi — 1.4 fm, xq = 0.6) are identical at 
the scale of the figure. As expected, a = 5 fm is too small (\Vn(o) /Vc(a)\ = 1.4). To 
obtain a good stability, radii larger than 6 fm should be used. 

The matching problem is illustrated in Fig. [51 where we show the wave function 
at E = 2 MeV, with a = 8 fm, N = 15. With the Lagrange functions, the matching 
between internal and external wave functions is quite smooth. For Sine functions, the 
matching is poor, which has a direct impact on the phase shifts. 

The accuracy of the different basis functions is illustrated in Table Q] where we 
compare various calculations of the phase shifts. The relative difference between left 
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Figure 5. 12 C+p I = wave functions with Lagrange and Sine functions at E = 2 
MeV (a = 8 fm, N = 15). Solid lines represent the internal wave function, and dotted 
lines the external parts. The Lagrange wave function is superimposed to the exact 
result. 



Table 1. 12 C+p I = phase shifts (in degrees) and matching parameters (|4.15|) for 
different bases (a = 8 fm). 



E (MeV) 




phase shift 




matching 


• parameter 


exact 


N = 7 


N = 10 


7V = 15 


N = 7 


N : 


= 10 


iV = 15 


Lagrange 


















0.5 


154.66 


112.90 


154.94 


154.59 


5.96 





65 


0.01 


1.0 


147.48 


144.22 


147.55 


147.48 


3.68 





28 


0.00 


1.5 


133.30 


311.02 


133.35 


133.30 


1.67 





20 


0.00 


2.0 


121.18 


299.30 


121.23 


121.18 


1.21 





16 


0.00 


Gaussian 


















0.5 


154.66 


179.36 


154.53 


154.53 


5.61 





02 


0.01 


1.0 


147.48 


146.52 


147.47 


147.47 


0.04 





01 


0.01 


1.5 


133.30 


130.28 


133.29 


133.29 


1.06 





00 


0.01 


2.0 


121.18 


112.15 


121.18 


121.18 


1.55 





01 


0.01 


Sine 


















0.5 


154.66 


2.03 


3.74 


6.68 


0.12 





66 


1.54 


1.0 


147.48 


28.90 


66.57 


109.98 


1.30 


4 


66 


29.35 


1.5 


133.30 


63.84 


94.53 


113.28 


11.46 


6 


97 


3.53 


2.0 


121.18 


71.76 


92.24 


105.34 


5.73 


3 


14 


2.38 



and right derivatives provides the matching parameter e defined as 
_ \^'(a)-u^'(a)\ 



«r t '(a) + < t '(a)]/2 
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This comparison is done for the different bases and at four typical energies. When e 
is small, the matching is obviously accurate. However, in some specific cases where 
M^ xt '(a) + wj nt '(a) is small, this parameter may be large, although the phase shift is 
fairly good. This parameter should therefore be considered as indicative only. Table [1] 
confirms that Lagrange and Gaussian functions represent a convenient basis, whereas 
Sine functions do not provide satisfactory results. Note that the use of Gaussian 
functions implies two additional parameters (bi, Xq) to define the basis. This choice has 
to be optimized for the different conditions and is therefore less direct than Lagrange 
functions. 

4-4- Application to a broad resonance: a + a 

The a + a system presents a well known rotational band based on the + ground state. 
The I = narrow resonance corresponds to the unstable ground state of 8 Be. In order 
to illustrate the R- matrix formalism applied to a broad resonance, we consider the I = 4 
partial wave. The experimental energy and width of the resonance are Er = 11.35±0.15 
MeV and T £s 3.5 MeV [97J. However, for broad resonances these properties may depend 
on the definition used (see section 3.6). A comparison between experiment and theory 
cannot be done with a high precision. 

The a + a potential of Buck et al. [98] has a simple Gaussian form, is /-independent, 
and reproduces the I — 0, 2, 4 experimental a + a phase shifts up to about 20 MeV. It 
is defined by 

V N (r) = -122.6225exp(-(r/2.132) 2 ), 

V c (r) = 4e 2 erf(r/1.33)/r. (4.16) 

Fig. [6] shows the phase shifts obtained with the Lagrange and Sine functions, with 
a = 8 fm. With Lagrange functions, the .R-matrix perfectly reproduces the exact phase 
shifts with N > 10. With Sine functions, the convergence is, as expected, much slower. 
Table |5] gives the phase shifts and e values under different conditions and at different 
energies. Again, the importance of the matching on the phase shifts is obvious. 

4-5. Application to a non-resonant system: a+ 3 He 

For I = 0, the experimental a+ 3 He phase shifts are well reproduced by the potential of 
Buck and Merchant [99], 

V N {r)= - 66.10 exp(-(r/2. 52) 2 ), 
V c (r) = 4e 2 (3 - {r/r c ) 2 )/2r c for r < r c , 

= 4e 2 /r for r > rc, (4-17) 

with rc = 3.095 fm. In Fig. [71 we use the Lagrange functions for two radii: a = 8 fm, 
and a = 5 fm. For a = 8 fm, the convergence is reached as soon as > 10. On the 
contrary the choice a = 5 fm is not consistent with one of the .R-matrix requirements: 




Figure 6. a + al=4 phase shifts with the Lagrange and Sine functions for different 
N values. For the Lagrange functions N > 10 is indistinguishable from the exact 
result. 



Table 2. a + a I = 4 phase shifts (in degrees) and matching parameters for Lagrange 
and Sine bases (a — 8 fm). 



E (MeV) 




phase shift 




matching 


• parameter 


exact 


N = 7 


N = 10 


N = 15 


N = 7 


N - 


= 10 


N= 15 


Lagrange 


















5 


0.79 


0.81 


0.80 


0.80 


0.05 





02 


0.00 


10 


24.46 


26.19 


24.58 


24.48 


0.42 





12 


0.00 


15 


125.88 


137.90 


126.00 


125.88 


0.31 





18 


0.01 


20 


141.55 


151.24 


141.53 


141.54 


1.26 


5 


67 


0.13 


Sine 


















5 


0.79 


0.53 


0.62 


0.68 


0.66 





65 


0.65 


10 


24.46 


15.55 


17.74 


19.84 


2.23 


2 


20 


2.17 


15 


125.88 


109.09 


116.76 


119.87 


2.16 


2 


09 


2.06 


20 


141.55 


139.82 


141.24 


141.37 


3.27 


5 


84 


6.21 



the nuclear contribution is not negligible. In that case, it is impossible to get a good 
convergence at all energies. 

In Fig. El we compare the wave functions at E = 8 MeV with the Sine and Lagrange 
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Figure 7. a+ 3 He i?-matrix phase shifts with Lagrange (a = 8 fm and a = 5 fin) and 
Sine (a = 8 fm) functions for I = 0. The thick lines represent the exact phase shifts. 

functions. In the former case, the matching between the internal and external wave 
functions is poor, which has an impact on the accuracy of the phase shifts. Conversely, 
the .R-matrix wave function with the Lagrange basis is indistinguishable from the exact 
wave function. 

Table [3] gives the phase shifts and matching coefficients under different conditions. 
For a = 8 fm, N > 10 provides reasonable phase shifts. If a increases (we choose here 
a = 10 fm), the .R-matrix approximation is of course still valid, but the number of basis 
functions must be increased. This is expected as the internal wave function must be 
reproduced in a wider range. In agreement with Fig. [51 the Sine functions provide a 
poor approximation of the phase shifts. 




Lagrange, a=8 fm 




Lagrange, a=5 fm 
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Figure 8. a+ 3 Hc I — wave functions with Sine functions at E = 8 McV (a = 8 
fm, N = 15). Solid lines represent the internal wave function, and the dotted line the 
external part with Sine functions. The Lagrange wave function is superimposed to the 
exact result. 



Table 3. a+ 3 Hc I = phase shifts (in degrees) and matching parameters for different 
bases. 



E (MeV) 








phase shift 




matching 


parameter 


exact 


N = 10 


N = 15 


N =10 


N = 15 


Lagrange 


a = 


8 f m 












5 






-63.97 


-65.02 


-63.96 


0.16 


0.01 


10 






-101.2 


-103.1 


-101.2 


0.89 


0.22 


15 






-125.9 


-128.9 


-125.9 


0.13 


0.01 


20 






-144.2 


-148.5 


-144.2 


0.14 


0.02 


Lagrange 


a = 


10 fm 












5 






-63.97 


-70.25 


-63.97 


0.30 


0.04 


10 






-101.2 


-109.5 


-101.3 


0.17 


0.02 


15 






-125.9 


-134.5 


-126.0 


5.19 


0.07 


20 






-144.2 


-149.4 


-144.3 


0.28 


0.02 


Sine, a = 


8 f m 














5 






-63.97 


-88.74 


-81.18 


2.66 


2.47 


10 






-101.2 


-127.0 


-117.2 


1.16 


0.81 


15 






-125.9 


-159.4 


-146.5 


2.66 


2.40 


20 






-144.2 


-183.7 


-165.4 


1.95 


1.82 



4-6. Application to a deep potential 

As mentioned before, the potential may contain additional (unphysical) bound states to 
simulate the Pauli principle. In the applications considered so far, the potentials involved 
a small number of forbidden states. The situation is different in heavy-ion reactions, 
where the number mi of forbidden states may be rather large. In such conditions the 
internal wave function presents several oscillations which must be reproduced by the 



The R-matrix Theory 



44 



basis functions to provide accurate phase shifts. 

This problem is illustrated here with the 13 C+ 12 C system (/ = 0, mi = 12). We take 
the deep Gaussian potential of |100j (depth of —273.7 MeV and range of 3.06 fm) for 
the nuclear interaction. The Coulomb potential has a point-sphere shape [see f!4.17p ]. 
with a radius of 8.2 fm. The i?-matrix phase shifts are determined with the Lagrange 
basis, and compared in Fig. [9] (upper panel) with the exact phase shifts. The channel 
radius is a = 10 fm. It is clear that small N values do not reproduce the phase shifts. 
Values larger than N ~ 35 are necessary to obtain accurate results. 

We show in Fig. [9] (lower panel) the wave function for N = 20,30,40 at E — 10 
MeV. Clearly, N = 20 is not able to give a good description of the oscillations of the 
internal wave function. The situation is of course improved with iV = 30, but only 
N = 40 is superimposed to the exact wave function. This problem is common to all 
potentials presenting many bound states. 

180 
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60 
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S -60 

to 
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Figure 9. Upper panel: 13 C+ 12 C I = phase shifts with the Lagrange functions 
(a = 10 fm). Lower panel: corresponding wave functions at E = 10 MeV. The exact 
wave function (superimposed on the N = 40 result) is shown as a thick line. 




E (MeV) 




4-7. Application to a non-local potential: e —H scattering in the static-exchange 
approximation 

The electron + hydrogen-atom scattering has been studied with various methods [5T] . 
In the simplest version, referred to as the static-exchange approximation |1011 1102] , the 
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problem is limited to a single-channel calculation, but involves a non-local potential. 
The local term U(r) reads 



We use here the traditional atomic units, involving the Bohr radius ao and the Hartree 
energy E = e 2 /a 

The non-local contribution is defined as 



where E is the total energy, S is the total spin of the electrons (S — or 1), and 
r > = max(r, r'). The l/r> singularity is responsible for numerical problems |1U3] , and 
deserves special attention in the literature. 

This example is used to illustrate the applicability of the Lagrange functions for 
non-local potentials (notice that the singlet phase shifts have been already determined 
in the i?-matrix framework [51]). Here the Gauss approximation associated with the 
Lagrange basis is expected to present convergence problems, owing to the singularity 
mentioned above. This problem will be addressed by performing, in parallel with the 
consistent Gauss approximation of the potential matrix elements [see (14.121) and (I4.13p ]. 
a "traditional" calculation using the Simpson method (with typically 1000 integration 
points) for the numerical quadrature. 

Results are shown in Tables 0] and for S = and S — 1, respectively, as a 
function of the wave number k. For the Lagrange-mesh method, we select a channel 
radius a = 14 ao, and consider various basis sizes. Comparing with the alternative 
method (Simpson integration) which is less sensitive to singularities in the potential, 
we immediately see that a reasonable accuracy (typically 0.02 rad.) can be achieved 
with small iV values at the Gauss approximation. However, when N increases the 
convergence is rather slow. If we do not use the Gauss approximation for the potential 
matrix elements, at least 5 digits are exact with iV = 20 and for a given a value. 
The sensitivity with respect to the channel radius is about 10 -4 rad. when going from 
a = 14 a to a = 16 ao- For comparison we give in Table H]the results of Burke et al. 
[5"T] . obtained for S = 0, up to k = 0.7 a^ 1 . 

Similar conclusions can be drawn for the S = 1 phase shifts, but the Gauss 
approximation is slightly more accurate because the antisymmetry of the spatial part 
of the wave function decreases the effect of the l/r> singularity. The 5 = and S — 1 
phase shifts are also in good agreement with those of Apagyi et al. |103] who use the 
Schwinger variational method |104j . much more sensitive to numerical problems than 
the i?-matrix method. 

4-8. Discussion of resonances 

In this section we determine the resonance energies and widths of the 12 C+p (/ = 0) and 
a + a (I — 4) systems, as an illustration of section EHfl In Fig. [TUl we show the functions 




(4.18) 




(4.19) 
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Table 4. i?-matrix phase shifts (in rad.) for e~— H scattering (5 = 0). The potential 
matrix elements are computed, either with the associated Gauss approximation (14.121) 
and (|4.13p . or with the Simpson quadrature. The results are compared with those of 
Burke et ai. [EE]. 









Gauss approximation, a = 14 


a 




N = 20, 


Simpson 




[SU 


k (giq 1 ) 


N 


= 20 


N = 40 


N = 60 


N = 80 


N 


= 100 


a = 14 a 


a = 16 ao 






0.1 


2 


3790 


2.3917 


2.3940 


2.3949 


2 


3952 


2.3959 


2.3958 


2 


3960 


0.2 


1 


8509 


1.8654 


1.8681 


1.8691 


1 


8695 


1.8703 


1.8702 


1 


8704 


0.3 


1 


4902 


1.5036 


1.5061 


1.5070 


1 


5074 


1.5081 


1.5081 


1 


5081 


0.4 


1 


2228 


1.2352 


1.2375 


1.2384 


1 


2388 


1.2394 


1.2395 


1 


2391 


0.5 


1 


0155 


1.0274 


1.0296 


1.0304 


1 


0308 


1.0314 


1.0315 


1 


0309 


0.6 





8532 


0.8650 


0.8673 


0.8681 





8684 


0.8691 


0.8690 





8685 


0.7 





7281 


0.7401 


0.7424 


0.7432 





7435 


0.7442 


0.7441 





7438 


0.8 





6347 


0.6470 


0.6494 


0.6502 





6506 


0.6513 


0.6512 






0.9 





5685 


0.5813 


0.5837 


0.5846 





5850 


0.5857 


0.5857 






1.0 





5251 


0.5384 


0.5409 


0.5417 





5422 


0.5429 


0.5429 







Table 5. See caption to Table S for 5 = 1. 









Gauss approximation 


, a = 14 


a 


N = 20, 


Simpson 


k (a, 1 ) 


N 


= 20 


N = 40 


N = 60 


N = 80 


N = 100 


a = 14 a 


a = 16 ao 


0.1 


2 


9080 


2.9077 


2.9077 


2.9077 


2.9077 


2.9075 


2.9076 


0.2 


2 


6799 


2.6794 


2.6793 


2.6792 


2.6792 


2.6791 


2.6792 


0.3 


2 


4622 


2.4614 


2.4612 


2.4612 


2.4611 


2.4611 


2.4611 


0.4 


2 


2588 


2.2576 


2.2574 


2.2574 


2.2573 


2.2572 


2.2573 


0.5 


2 


0721 


2.0706 


2.0703 


2.0702 


2.0702 


2.0701 


2.0701 


0.6 


1 


9032 


1.9013 


1.9009 


1.9008 


1.9007 


1.9006 


1.9006 


0.7 


1 


7521 


1.7497 


1.7492 


1.7490 


1.7490 


1.7488 


1.7488 


0.8 


1 


6180 


1.6150 


1.6145 


1.6143 


1.6142 


1.6140 


1.6141 


0.9 


1 


4998 


1.4963 


1.4957 


1.4954 


1.4953 


1.4951 


1.4952 


1.0 


1 


3959 


1.3919 


1.3911 


1.3909 


1.3908 


1.3905 


1.3905 



R(E) and 1/S(E) for both systems. The resonances energies are calculated from the 
roots of equation (I3.42p . i.e. from the crossing points of both curves. The numerical 
values, as well as the widths (I3.45[) are given in Table [6] 

The 12 C+p (/ = 0) state is a typical example of a narrow resonance. The phase 
shifts near 0.42 MeV are well described by the Breit-Wigner approximation. Although 
the R matrices are rather different for a = 8 fm and a = 10 fm (see Fig. [TO]) , the energy 
and width are weakly sensitive to the channel radius (see Table [6]). This result is typical 
of narrow and isolated resonances. 

The situation is different for the a + a (1 = 4) system, which is typical of a broad 
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Figure 10. i?-matrix Ri{E) and inverse of the shift function Si(E) for the 12 C+p 
(I = 0) and a + a (1 = 4) collisions. Solid lines correspond to a = 8 fm and dotted 
lines to a = 10 fm. The circles show the roots of (|3.42|) . 



resonance (around 12.5 MeV, see Fig. [6j). The resonance properties are rather sensitive 
to the channel radius, whereas the phase shifts are almost identical. Table |U] shows that 
the second eigenvalue might be related to the 4 + resonance. However the width is not 
significantly smaller than the energy difference with neighbouring states. In that case, 
a single-level approximation is not well adapted to reproduce the phase shifts in a wide 
energy range. 

Table 6. Resonance energies and widths (in MeV) of the 12 C+p (I = 0) and 
a + a (1 = 4) systems. The Lagrange basis is used with N = 15. 





a = 


8 f m 


a = 


10 fm 


Eigenvalue 


Er 


r 


Er 


r 


12 C+p (/ = 0) 










1 


0.418 


0.0375 


0.417 


0.0369 


a + a (I = 4) 










1 


6.92 


5.14 


4.67 


3.68 


2 


12.3 


2.47 


10.2 


2.74 


3 


19.8 


7.08 


14.0 


3.33 
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4-9. Application to bound states 

The R matrix is applied to the 13 N ground state (J w = l/2~,l = 1) described by a 
12 C+p potential model. In order to simulate the experimental binding energy (—1.94 
MeV), the depth of the Gaussian potential used in section I4T31 for I = has been modified 
to -55.3 MeV. 

The calculation is performed with the Lagrange basis, and the results are presented 
in Table [7J With a = 8 fm, N = 7 does not provide the exact energy. This situation 
is of course improved by increasing N. For N = 15, we have 5 exact digits, and the 
matching parameter is close to 0. In each case, 4 iterations are sufficient in system (13.541) 
to provide an energy stable by better than 0.1 keV. The ANC value Ci [see ( 13. 52ft ] is 
also very stable. 

As in section 4.3 for the phase shifts, the choice a = 5 fm is too small to make the 
nuclear interaction negligible. Although the matching parameter is quite acceptable, the 
convergence is slower and the final result does not converge to the exact energy. This 
means that a small matching parameter is not sufficient to ensure that the calculation 
is accurate. Testing the stability against the channel radius is a more severe test. 

Table 7. Energy of 13 N ground state (in MeV) with a Lagrange basis and different a 
and N values. The exact binding energy and ANC are —1.942 MeV and 2.063 fm -1 / 2 , 
respectively. 

Iteration a = 8, iV = 7 a = 8,iV = 10 a = 8, iV = 15 a = 5,iV=10 



1 


-2.012 


-2.052 


-2.053 


-3.236 


2 


-1.897 


-1.940 


-1.941 


-1.763 


3 


-1.899 


-1.941 


-1.942 


-1.894 


i 


-1.899 


-1.941 


-1.942 


-1.881 


5 


-1.899 


-1.941 


-1.942 


-1.882 


6 








-1.882 


(fm^ 1 / 2 ) 


2.060 


2.071 


2.072 


2.034 


e 


9.1 


0.19 


0.005 


0.002 



4.. 10. Application to a multichannel problem: a+d 

Here we use the «+deuteron potential of Dubovichenko |105] to investigate the R- 
matrix formalism in a multichannel problem. This nuclear potential Vjv(r) contains 
central V c (r) and tensor V t {r) forces, defined as 

V N (r) = V c (r) + V t (r)S 12 

yjj) =Voexp(— ar 2 ) 

V t {r) =V 1 exp(-(3r 2 ) J (4.20) 

where 

S 12 = ^(S-r) 2 -2S 2 (4.21) 
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is the usual tensor operator. The Coulomb potential is the bare potential. With 
V = -91.979 MeV, V x = -25.0 MeV, a = 0.2 fin" 2 , = 1.12 fm" 2 , this potential 
reproduces most of the 6 Li ground-state properties, in particular the binding energy, 
the quadrupole moment, and the c?-wave admixture amplitude. The spin and parity of 
the deuteron being 1 + , non-natural parity states (i.e. it = (— 1) J+1 ) involve two I values: 
I — | J — 1| and / = J+1. The system to be solved (13.77P involves the potentials (we 
use the notation of section I3.10p 

Vn = V C + V C - 2( J - l)V t /(2J + 1) 



V 12 = V 2 i = 6y J(J + l)V t /(2J + 1) 

V 22 = V c + V c -2(J + 2)V t /(2J+l) (4.22) 

As in previous sections the .R-matrix phase shifts are computed in various 
conditions, and compared to the "exact" results obtained by the Numerov algorithm. 
The .R-matrix calculations are performed with Lagrange functions. To compare with 
experiment (106J, the collision matrix is diagonalized [see (j2.37|) ] and parametrized as 

u= ( cose -sine \ t exp(2z<5i) \ / cose sine \ 
\ sine cose J \ exp(2i5 2 ) J \ — sine cose J ' 

where Si and 82 are the (real) eigenphases, and e is the mixing angle (removing the 
tensor force provides e = 0). These values are presented in Fig. [TT] for J n = 1 + and 
compared with the exact solutions (also shown in Figs. 4 and 5 of |105j ). We choose 
a = 8 fm and show the result for N = 7. With this small number of basis functions, 
slight deviations can be observed. As soon as N > 7, the .R-matrix phase shifts cannot 
be distinguished from the exact values. 

Other applications of the .R-matrix theory to multichannel problems can be found 

in 



4-11. Application to propagation methods 

To illustrate the propagation of the R matrix, we present here a numerical example, 
with the a + a system. We choose a = 80 fm, which is of course unnecessarily large for 
this system, but typical of long range potentials, where propagation methods should be 
used. The goal of this example is just to provide a numerical illustration of the method. 
The R matrix at 80 fm has been computed in 3 ways: (i) exactly with the Numerov 
algorithm, (ii) with 100 basis functions without propagation (N$ = 1), (hi) with 100 
basis functions split in several intervals. The .R-matrix basis functions are the Lagrange 
functions, but any other choice would provide similar conclusions. Matrix elements over 
an interval (ai,a 2 ) are given in Appendix C [see (!C14p - (lU18p . 

Table [8] gives the I = .R-matrix at a = 80 fm for typical energies. Without 
propagation (N$ = 1), we reproduce at least 4 significant digits. This is still true for 
Ns = 2, but the computer time is reduced by more than a factor of two. Essentially the 
difference is that, in the former case we have to invert one 100 x 100 matrix, whereas 
the latter calculation requires two inversions of 50 x 50 matrices. The computer time is 
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Figure 11. a+d cigenphascs (upper panel) and mixing angles (lower panel) for 
J T = 1 + . The solid lines represent the exact calculation and the dotted lines 
correspond to the i?-matrix calculation with a — 8 fm, N = 7. The experimental 
data are taken from |106j . 



still reduced with N$ = 4 (four inversions of 25 x 25 matrices). The precision is however 
reduced since the first interval (0,20 fm) only contains 25 basis functions. Increasing 
this number to 35 (and hence the full basis size to 110 functions) provides the exact 
i?-matrix (with at least 4 significant digits) with similar computer times. 



Table 8. i?-matrix for the a + a collision (I = 0) and for different N$ values 
(N = 100, a = 80 fm). For each interval the number of basis functions is 100 /Ns- 

E (MeV) exact N s = 1 N s = 2 N s = 4 

IT 2.648 x 10~ 2 2.648 x 10" 2 2.648 x 10~ 2 2.643 x 10~ 2 

10 -3.729 x 10" 1 -3.729 x 10" 1 -3.730 x lO" 1 -3.957 x 10" 1 

15 1.195 x 10~ 3 1.195 x 10~ 3 1.195 x 10~ 3 1.157 x 10~ 3 

20 -1.246 x 10~ 2 -1.246 x 10~ 2 -1.246 x 10~ 2 -1.262 x 10~ 2 

time (ms) 1.2 0.50 0.25 



4-12. Application to capture reactions: 12 C(p,~f) 13 N 

The 12 C(p,7) 13 N reaction is the first reaction of the CNO cycle. It represents an ideal 
candidate for an .R-matrix treatment since the l/2 + (/ = 0) resonance at Er = 0.42 
MeV determines the 5-factor ( 12. 55ft in a wide energy range. As in previous sections, our 
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goal here is not to find the best fit to the data (see p2]). Rather, we want to illustrate 
the different contributions to the capture matrix elements ( I3.6ip . and to discuss various 
approximations. Our procedure is as follows: 

(i) The initial states are determined as in section 14.31 

(ii) As our goal is to illustrate the different contributions to the capture matrix elements, 
the bound-state wave function (and the corresponding ANC) is computed exactly 
with the Numerov algorithm [90]. This avoids the optimization of the basis for the 
bound state. 

(iii) The basis functions for the initial state are chosen as Lagrange functions (N = 15), 
and we determine the internal and external contributions from (I3.62p and ( I3.63p . 
respectively. 

The matrix elements are computed as in section |3T8| but the i?-matrix expansion of 
the final state is replaced by the exact wave function. Table Ogives, at typical energies, 
the exact values of the matrix elements (i.e. with scattering wave functions obtained 
from the Numerov method), as well as their values in the i?-matrix theory with a = 8 
fm and a = 10 fm. Large values must be used to ensure that the nuclear interaction is 
negligible. Several comments can be made: 

(i) As expected the matrix element presents a maximum at E — 0.42 MeV, as the 
initial potential has been fitted to provide a resonance at this energy. 

(ii) Since the external scattering wave function (I3.63f) involves the phase shift, the 
external contribution also presents a resonant behaviour. 

(iii) Each term in (13.61)) depends on the channel radius. Their sum, however, should be 
insensitive to its choice, if the conditions of the calculation are properly defined. 
The relative difference is maximum near the resonance but is always less than 1%. 

For the sake of completeness we show in Fig. [T2]the 12 C(p, 7) 13 N S'-factor compared 
to the available data sets [1071 1108) . As usual in the potential model, a spectroscopic 
factor S should be introduced for a more realistic comparison with the data. The 
spectroscopic factor scales the total cross section. Values lower than unity mean that 
the final ground state is more complicated than a simple 12 C+p structure. With 5=1 
the theoretical 5-factor is larger than the data. A reasonable agreement can be obtained 
with S = 0.45, although the high energy part is slightly overestimated by the model. 
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Table 9. 12 C(p, 7) 13 N matrix elements [see (g^H and ([535]) ] with La grange functions 
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2002.1 
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2412.7 
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experimental data are from [107] (open circles) and |108| (closed circles). 
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5. The phenomenological R matrix 

5.1. Introduction 

The goal of the phenomenological .R-matrix method is to use a parametrization based 
on expression (I3.23P of the R matrix or its multichannel generalization (j3.86|) with a 
finite number of poles. The properties of these poles are adjusted to some data, in place 
of being derived from some Hamiltonian, as in the calculable approach. We present here 
various applications in nuclear physics. In particular, this technique is very successful 
in nuclear astrophysics |109j . where the main issue is to fit cross-section data, and to 
extrapolate them down to stellar energies at which direct measurements are in general 
impossible (see, for example, [21]). Another recent application is the analysis of low- 
energy scattering data in experiments involving radioactive beams. Those experiments 
usually probe the nuclear structure at low level density, well adapted to the .R-matrix 
formalism (see, for example references |11U1 1112] for 8 B+p, 11 C+p and 18 Ne+p, 
respectively). The method is of course not limited to elastic scattering, but can be 
extended to inelastic |113j as well as to transfer |114] reactions. 

One of the main drawbacks of the phenomenological .R-matrix formalism is that, 
though the pole energies and reduced widths are associated with physical properties, 
they cannot be directly compared with experiment. In particular, the .R-matrix 
parameters significantly depend on the channel radius. This is in contrast with the 
so-called "observed" data, such as resonance energies and widths directly fitted to 
experimental data which are by definition independent of such a radius. However, these 
"observed" values of the energy and width of a resonance depend on the assumptions 
made about the theoretical description of the resonance (see section for a part of the 
possible definitions). Within the phenomenological R- matrix approach, the "observed" 
quantities have a rather precise definition from which the .R-matrix parameters can be 
determined. For an isolated resonance in some partial wave /, the observed energy Er 
and width Tr are obtained in the single-channel approach by fitting the elastic cross 
section (12.461) where the phase shift is parametrized by the generalized Breit-Wigner 
expression ( 13.480 . 

-T 

5f w + arctan 2 * . (5.1) 
Er — E 

This expression may be simplified by neglecting <pi or made more realistic by multiplying 
Tr by Pi(E)/Pi(Er). Notice that, except for the simplest Breit-Wigner approximation, 
these expressions slightly depend on a through <pi or The "observed" reduced width 
7^ bs is then extracted from Y R with the relation similar to (13.500 . 

T R = 2 1 l bs P l {E R ). (5.2) 

The corresponding "formal" parameters, i.e. the pole location E n i and the reduced 
width 7^ can then be deduced for a given channel radius a but their determination is 
not immediate because of the shift factor Si and its energy dependence. The problem is 
more complicated if several resonances or several channels must be taken into account. 
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In the next sections, we present different methods to link the "formal" .R-matrix 
parameters with the "observed" values. The simplest and most common case deals 
with a single isolated resonance, described by one pole (section 15.21) . However some 
applications require to include several poles in the same partial wave or, in other words, 
to consider interference effects between these poles (section 15. 3p . Other applications 
require to take several channels into account (section I5.4p . As a first application, 
we consider the 12 C+p elastic scattering from 0.3 MeV to 1.8 MeV, where accurate 
data exist for many years (section 15. 51) . Then more recent multichannel application 
to 1 H( 18 Ne,p) 18 Ne(g.s.) and 1 H( 18 Ne,p') 18 Ne*(2+, 1.887 MeV) will be shown as an 
example of radioactive-beam experiments (section 15. 6f) . Finally, the .R-matrix method 
will be applied to radiative capture reactions, with 12 C(p, 7) 13 N as a typical example 
(section 15.71) . 



5.2. Single-pole approximation of elastic scattering 

In section 13.61 we have presented a general procedure to define the resonance energy 
and width from the .R-matrix expression. Let us particularize (I3.23P to a single pole 
with energy E\ and reduced width 7 2 (we drop the index I for the sake of clarity); the 
R matrix is then given by 

Ri(E) = z^z- (5-3) 

This approximation is frequently used at low energy, where single isolated resonances 
are present. The .R-matrix phase shift associated with ( 15.31) is given by (I3.47p . To go 
further, let us use the Thomas approximation [S], which consists in a linearization of 
the shift function Si(E). Near the pole energy E±, it reads 

Si(E) « S^E,) + (E- E 1 )S' l (E 1 ), (5.4) 

where is the derivative of the shift factor with respect to energy (which appears 
both in the wave number k and in the Sommerfeld parameter rj). The validity of this 
approximation is supported by Figs. [TJ [2]and[31 where it is clear that, in a limited energy 
range, the linearization of the shift function is quite appropriate. 

Comparing expression ( I3.47P of the phase shift at the Thomas approximation with 
the Breit-Wigner expression (15. ip . one derives the observed properties (En, Jobs) from 
the formal parameters (£1,71) of an isolated pole. The observed energy Er reads 

= < 5 ' 5 ' 

where the shift between both energies depends on E\. This shift also depends on a. It 
is in general non-negligible, unless 7 2 is very small. The observed reduced width reads 

l2 ° hs = TT^spo ' (5 ' 6) 

It also depends on E±. 
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In practice, however, the reversed relationships are needed. Indeed, in many cases, 
observed values (Er, 'jobs) are known from experiments (or from other theoretical works) 
and one wants to derive the corresponding i?-matrix parameters (#1,71). The inverses 
of (I5.5p and (15.61) are obtained by linearizing the shift factor Si(E) around Er as 

E 1 = E R + >HSi(E R ). (5.8) 

These results are well known [5]. 

In the literature, the reduced width is frequently expressed in units of the Wigner 
limit [BH] 

* = £• (5 - 9) 

which provides the dimensionless reduced widths 

^ = 7! 2 /7^- (5.10) 

Notice that Lane and Thomas [5] define the dimensionless reduced width as 9 2 = 
j 2 bs fj,a 2 /h 2 without an explicit reference to the Wigner limit. In general, a value of 
9 2 not far from unity indicates the occurrence of a cluster structure, i.e. the colliding 
nuclei partly conserve their identity within the resonance. Conversely, a plausible guess 
on 9 2 (possibly inspired by a model) can be used to estimate the width of a resonance. 
This concept of dimensionless reduced width, as defined by (15.1 Oft , was first used by 
various authors |115l 11161 1117] in the interpretation of low-energy scattering data. 

The difference between observed and formal parameters is illustrated by numerical 
applications in Table fTUl for the narrow l/2 + resonance in 12 C+p (Er = 0.42 MeV, 
T R = 32 keV) and for the broad 1" resonance in l2 C+a (E R = 2.42 MeV, T R = 420 
keV). The formal parameters (£1,71) are calculated for various channel radii and are 
then used to determine the phase shifts shown in Fig. [13] As usual in applications of 
the phenomenological i?-matrix method, the channel radii are somewhat smaller than 
in the "calculable" variant. 

From Table EH it is clear that the dimensionless reduced widths 9 2 bs are less 
dependent on the channel radius than / y 2 bs . For both systems, and in particular for 
12 C+a, the 9 2 bs values are rather large. These states can therefore be considered as 
cluster states. As expected, the formal parameters ( 15. 7ft are strongly dependent on the 
channel radius. However, the corresponding phase shifts in Fig. [13] are very close to 
each other, in particular for the narrow resonance in 12 C+p. For the broad resonance in 
12 C+a, the phase shift is rather stable near the resonance energy, but more significant 
differences appear at higher energies. In such a case, the validity of the Breit- Wigner 
approximation is more limited. In other words, the single-pole approximation (15.31) 
should be replaced by an R matrix containing several terms. 
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Table 10. i?-matrix parameters (|5.7[) and (|5.8p for resonances in 12 C+p and 12 C+a 
(in MeV). The observed reduced widths are obtained from (|5.2j) . 



C+p (J* = 1/2+, Z = 0, E R = 0.42 MeV, T R = 32 keV) 
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Figure 13. 12 C+p (I = 0) and 12 C+a (I = 1) phase shifts in the single-pole i?-matrix 
approximation computed with the parameters of Table ITDl for different channel radii. 



5.3. Multiresonance elastic scattering 

As mentioned above, the single-pole approximation is often valid at low energies. 
However, for nuclei with a high level density, several resonances may appear with 
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the same angular momentum and parity. Then the interferences between different 
resonances and/or bound states in the same partial wave may be important. Typical 
examples are 12 C+a or 14 N+p where several partial waves involve more than one state, 
even in a limited energy range. 

In this case, the link between formal and observed parameters is more complicated, 
and was addressed for many years in a rather indirect way (see, for example, [HI 1118] ). 
More recently, an iterative procedure was proposed to determine the formal i?-matrix 
parameters from observed values [23] • This method is valid for single- channel systems 
only. A generalization to multichannel problems was then developed by Brune |119j . 
The idea is to propose an alternative parametrization of the R matrix, where the input 
parameters are the observed data. It is based on the invariance of these values when 
the boundary parameters B c [see (I3.80p ] are changed. We briefly summarize Brune's 
method here, by assuming elastic scattering only. A more general presentation can be 
found in [119] . 

Let us assume N resonances in a given partial wave, with observed energies Em 
and widths T^. For each resonance, an observed reduced width is defined, according to 
052]), as 

7 L,, = rV2Pi(£*i). (s.ii) 

In the notations of |119j . ( 15.71) is written as 

1 - lobs,i b A^Ri) 

From these expressions, the formal pole energies E n , used in the iV-pole i?-matrix 
expansion, are obtained from 

Nb n = E n Mb n (5.13) 

where matrix elements iVy and My are given by 

m \ Em + 7fSi(Em) for i = j, 

[ 7j 13 E m -E R] iul ^ T h 

and 

{1 for % — j, 

~ ~ MEnJ-StlEy) { , (5.15) 

The generalized eigenvalue problem (15.131) provides the formal energies E n . The formal 
reduced-width amplitudes 7 n are derived from the eigenvectors b n as 

N 

ln = Y, h njl3- (5-16) 

i=i 

When = 1, it is easy to see that ( 15.7ft and ( 15. 8j) are recovered. This method provides 
an efficient way to derive the i?-matrix parameters. It represents the starting point of 
an alternative parametrization of the R matrix, proposed by Brune |119] . 
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This formalism is illustrated in Table [TT] and Fig. [TH with the 14 N+p system. In the 
J w = 3/2 + partial wave (1 = 0), we take three resonances into account. The observed 
values Em and T R i [91] are given in Table [TTJ As for the single-pole approximation, 
the formal parameters (E n ,~f n ) do depend on the channel radius. Fig. [T41 shows that 
the influence of a on the corresponding phase shifts is weak near the narrow resonance 
at 0.987 MeV (this resonance presents a small 6 2 value and has thus a complicated 
structure). As expected, it is more important in the vicinity of the broader states at 2.2 
and 3.2 MeV. 

Table 11. i?-matrix parameters (in MeV) with 3 poles for the J x = 3/2 + ,Z = 
partial wave in 14 N+p. 



observed values 




a = 5 fm 


a = 6 fm 


a = 7 fm 


E m = 0.987 


lobs,l 


0.0082 


0.0054 


0.0039 


T R1 = 0.00367 


°obs,l 


0.0031 


0.0029 


0.0028 






0.0097 


0.0061 


0.0042 




E 1 


0.981 


0.983 


0.985 


E m = 2.187 


lobsfl 


0.114 


0.086 


0.069 


T R2 = 0.2 


°obs,2 


0.043 


0.047 


0.050 




ll 


0.119 


0.089 


0.070 




E 2 


2.153 


2.165 


2.172 


E R3 = 3.209 


lobs, 3 


0.053 


0.042 


0.034 


T R3 = 0.14 


ft" 2 

u obs,3 


0.020 


0.023 


0.025 




ll 


0.051 


0.041 


0.034 




Ez 


3.199 


3.202 


3.204 




Figure 14. 14 N+p i?-matrix phase shifts (J n = 
of Table [IU 



3/2 + ) computed with the parameters 
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5-4- Phenomenological parametrization of multichannel collisions 

To deal with inelastic or transfer cross sections, the .R-matrix formalism must be 
extended to several channels. Rigorously this is even true in elastic scattering involving 
particles with non-zero spins. When the spins are different from zero, several (II) values 
[see (12.29jl ] are to be considered, each of them being characterized by a partial reduced 
width. 

Let us consider the multichannel R matrix ( I3.86p . In the calculable variant, the 
number of poles is equal to the product of the number of basis functions N by the 
number of channels. In the phenomenological approach, including all partial waves is 
often not realistic, since too many parameters may be involved. At low energies, the 
choice of the relative angular momenta / is guided by the penetration factor, i.e. the 
minimum / value for given J and ir is often adopted. Such a selection mechanism does 
not exist for the channel spin I. In general, all possible values should be considered. 
The choice of the relevant partial waves is made according to the quality of the fit. 
For these reasons, the multichannel approach in the phenomenological variant of the R 
matrix is most often restricted to two channels. 

Here we limit the presentation to a single pole at energy Ei, the most frequently 
used approximation. In that case, the R matrix (I3.86P is reduced to 

RcAE) = (5.17) 

where 7 C is the reduced-width amplitude of the pole in channel c (for simplicity we drop 
the index 1 of the reduced widths). These parameters are real. In this approximation, 
the property 

R 2 CC ,(E) = R CC (E)R C , C ,(E), (5.18) 
provides the multichannel Z matrix (13.891) as 



(ak r ,)- ll2 O r , 



Sec 1 — L c i \ R cc R c i c i 



With the help of (IB2[) . we obtain the inverse as 

L c i yj R CC R C 



(ak c ) l ' 2 0~ l 



0~cc' + 



(5.19) 



(5.20) 



1 — RccL c 

where the logarithmic derivatives L c (I3.28j) and outgoing functions O c are defined for 
each channel. Using definition (I3.88P for a real R matrix, one obtains 

V P c Pc'RccRcU 



2i- 



(5.21; 



1 — J2c RccL c 

where <p c is the hard-sphere phase shift (I3.26P in channel c. 

From the denominator in (15.2 ip . the resonance energy and the reduced width are 
direct generalizations of (I5.5P and ( 15.61) . 

TP - TP ^2dc S c(El) , s 

E R - 2QKT? V ( b - 22 > 

l + EcX S c( E l) 



lobsfi— -I . ^ ^2QI(T?\' (5.23) 
l + Eclc^c(El) 



2 
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where index I in the penetration and shift factors is replaced by the more general channel 
index c. With (15.17p . the collision matrix is parametrized at the generalized Breit- 
Wigner approximation [96] as 



0~cc> + 



iJT c (E)T c ,(E) 



(5.24) 



Er — E — iT{E)/2 

where r c is the observed partial width in channel c and V = J2 C T c - A simple calculation 



gives 



T e (E) = 2^ obStC P c (E) 



(5.25) 



and r c ^ = T c (Er). Again, the .R-matrix parameters (Ei,%) can be deduced from the 
observed values as 



7c 



7o6s,c 



E 1 = E R + Y.'TiSc{E R ) 



(5.26) 



(5.27) 



These analytical formulas are direct extensions of ( 15. 7ft and ( 15. 8ft obtained for elastic 
scattering in the single-pole approximation. 



5.5. Application to the l2 C+p elastic scattering 

Data on 12 C+p exist for many years. In particular, elastic scattering cross sections 
have been measured with a high precision [120J. This system is well adapted to the 
phenomenological .R-matrix approach, since the level density near threshold is quite 
low. This example should be considered as a typical application of the method. Similar 
fits have been done on the 16 0+p I = elastic phase shifts |121j and on the 1 O+p cross 
section [48J. 

From [S3], three resonances are expected in the energy range covered by the data: 
1/2+ at 0.421 MeV, 3/2~ at 1.558 MeV and 5/2+ at 1.603 MeV. The pole corresponding 
to the bound state is neglected. Data sets at three cm. angles {6 = 89.1°, 118.7°, 146.9°) 
are available. The smallest and largest angles are fitted simultaneously by using the 
single-pole approximation fl 5 . 3 [) for the resonant partial waves. For other partial waves, 
the hard-sphere phase shift is used. This is consistent with the absence of resonance 
(Ri = 0), but plays a minor role in the cross sections. Replacing the hard-sphere phase 
shifts by zero provides essentially the same fits. 

The observed resonance properties are given in Table fT2l for different channel radii. 
Clearly the results are almost independent of a, as expected from physical arguments. 
The fitted values are consistent with the literature |120j . and the corresponding cross 
sections are shown in Fig. [15] for both scattering angles. The three channel radii provide 
fits which are indistinguishable at the scale of the figure. As expected |120j . the R- 
matrix parametrization reproduces the data very well, not only in the vicinity of the 
resonances, but also between them, where the process is mostly non-resonant. This 
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technique is very successful in the analysis of recent data involving radioactive beams 
(see, for example, [TTUl ITTT1 ITT2] ). 



Table 12. i?-matrix parameters from a simultaneous fit of 12 C+p scattering data |120j 
at = 89. 1° and 146.9°. Resonance energies Er are expressed in MeV and widths Tr 
in keV. 

7^ = 1/2+ J^ = 3/2" J w = 5/2 + 
Er T/j E r T r E r T r 
0.427 33.8 1.560 51.4 1.603 48.1 
0.427 32.9 1.559 51.4 1.604 48.1 
0.427 30.9 1.558 51.3 1.606 47.8 
0.424 33 1.558 55 1.604 50 



a = 4 fm 
a = 5 fm 
a = 6 fm 
Exp. [120] 




5.6. Application to the 18 Ne(p,p') 18 Ne(2 + ) inelastic scattering 

We present here an application of the phenomenological i?-matrix method to inelastic 
scattering. The 18 Ne(p,p') 18 Ne(2 + , 1.887 MeV) cross section has been measured in 
parallel with elastic cross sections |113j . These data were obtained at various angles and 
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complemented a previous data set, obtained at lower energies, and aimed at investigating 
elastic scattering only |122j . 

Here our goal is not to repeat the analysis of |113j . where several angles were 
simultaneously included, but where previous elastic data were not considered. Instead, 
we select a single angle but cover a broader energy range by including data sets of 
|122[I113] in a global fit. Both experiments measured the elastic cross sections in different 
energy ranges, but also at slightly different angles. We select the elastic data sets of 
|113] at Oiab = 6.2° and of |122] at 9i a b = 4.9°. As the angular dependence is weak, we 
combine these both data sets at a common angle, taken as the average {9i a b — 5.6°). 
For the inelastic cross section, the experimental angle [1 13j 6i a b = 6.2° is used. 




a 
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10 
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() () 
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Figure 16. i?-matrix fits of 18 Ne(p,p) 18 Ne elastic (upper panel) and 
18 Ne(p,p') 18 Ne(2 + ) inelastic (lower panel) cross sections with the parameters of Table 
IT51 The data are from [122] (full circles) and |113| (open circles). The fits are done 
with a = 4.5 fm (dashed lines), a = 5.0 fm (solid lines), and a = 5.5 fm (dotted lines). 



The experimental data and the corresponding i?-matrix fits are presented in Fig. [TBI 
As suggested in |113j . the fits are performed by including three resonances (J n = 
l/2 + , 5/2 + , 3/2 + ), which are characterized by their energy Er and their partial widths 
Ti and r 2 corresponding to the p+ 18 Ne(0 + ) and p+ 18 Ne(2 + ) channels, respectively. 
The fitted parameters are given in Table [13] for different channel radii. As expected 
we essentially reproduce the results of |122[ 1113] . The l/2 + resonance at 1.06 MeV is 
below the inelastic threshold (r 2 = 0), and corresponds to a single-particle state, with 
a large reduced width. The higher-lying resonances (5/2 + ,3/2 + ) correspond to s states 
in the p+ 18 Ne(2 + ) channel. They present a dominant width in that channel (r 2 3> fi) 
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and correspond to a significant fraction of the Wigner limit ( 15.9ft . These resonances 
are hardly visible in the elastic data, and could not be properly analyzed without the 
inelastic cross sections. 

Owing to the use of a radioactive beam, the error bars are rather large in the 
inelastic cross sections, and the sensitivity of the .R-matrix parameters to the channel 
radius is slightly stronger than in 12 C+p. In .R-matrix analyses, this sensitivity should 
be taken into account in the evaluation of the recommended error bars. 

Table 13. i?-matrix parameters from a simultaneous fit of elastic and inelastic 18 Ne+p 
cross sections (see text). Resonance energies are given in MeV and widths in keV. The 
bracketed values represent the dimensionless reduced widths. 







a = 4.5 fm 


a = 5 fm 


a = 5.5 fm 


1/2+ 


Er 


1.064 


1.063 


1.062 




Ti 


101 (0.26) 


98 (0.23) 


94 (0.20) 


5/2+ 


Er 


2.773 


2.768 


2.766 




Ti 


7 (0.01) 


6 (0.01) 


7 (0.01) 




r 2 


85 (0.46) 


85 (0.40) 


80 (0.34) 


3/2+ 


Er 


3.108 


3.080 


3.036 






11 (0.01) 


11 (0.01) 


9 (0.01) 




r 2 


300 (0.47) 


379 (0.57) 


306 (0.48) 


x 2 




0.43 


0.44 


0.54 



5.7. Radiative capture reactions 

5.7.1. Extension of R-matrix formalism 

The general formalism of radiative-capture reactions has been given in section 12.41 
In the .R-matrix theory, the matrix element U ( I2.58P is split in internal and external 
contributions as [21] 



U = U int + U ext , (5.28) 

where we have dropped all indices for the sake of clarity. The terms U int and U ext of the 
r.h.s. involve the internal and external wave functions, respectively. We assume here 
single-channel calculations, but the spins may be different from zero. 

By using expansion ( 13. 36ft for the initial radial wave function uf^ and = 1, the 
internal part of (12.581) becomes 



fr e ^" f) v ^(E)r n (E) 

int ~\l-L k R^\Y n E n — E ' {b - ZJ} 
where we have introduced the formal gamma width of pole n as 



r 



7,n 



2J f + 1 8tt(A + l)^ A n +1 
2J + 1 A(2A + 1)!! 2 



(V J f*f\\M aX \\$Z*)~ (5.30) 



2 
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with k 1>n = (E n — Ef)/hc. The matrix element is calculated over the internal region 
only. Its definition involves a state $^ 7r corresponding to the pole E n , whose radial part 
r _1 0^ 7r (r) is defined like in the orthogonal basis (13. 22ft but for partial wave Jn. Matrix 
element (15.291) depends on the particle width T n (E) (13.451) and on the gamma width at 
collision energy E, 



(for simplicity, V ltn {E n ) is often denoted as T 7jn ). In (15.291) . e n = ±1 is the product 
of the signs of the matrix element in ( I5.30P and of the reduced width amplitude. For 
a single-pole approximation, this sign does not play a role but determines interference 
effects in multi-pole calculations. 

The external contribution is determined as in (13.631) . A calculation similar to the 
one leading to (I2.64p provides 



E. This contribution is often referred to as "direct capture". In fact it is closely related 
to the internal term through the collision matrix. A resonant behaviour of the collision 
matrix affects the external term (15.321) (see also Table [H]). As usual, in the calculable R- 
matrix theory, the total matrix element (I5.28f) should not depend on the channel radius 
a, although each individual term does depend on a. In addition, one easily shows that 
both terms present an identical phase factor. The calculation can therefore be reduced 
to real expressions. 

In the calculable i?-matrix, the gamma width (I5.30p is computed from basis 
functions. In the phenomenological variant, the constant Y 1)n (E n ) appearing in (15.311) 
becomes a parameter. The treatment of radiative-capture reactions therefore requires 
one additional parameter for each pole, the gamma width r i}n (E n ) (and the associated 
interference sign), and a global parameter, the ANC C^ f7Tf of the final bound state. The 
latter parameter is sometimes available independently. As for elastic widths, the fitted 
values of the gamma widths may depend on the channel radius. The importance of this 
dependence will vary with the amplitude of the external contribution. 

5.7.2. Isolated resonance approximation 

Let us consider the single-pole approximation (15. 3p . Starting from expression (15.291) 
of the internal matrix element U- intl a simple calculation using definition ( I3.34p of the 
phase shift provides near the resonance energy Er the approximation 




(5.31) 




where C l /7r/ is the ANC of the final bound state and Uf. n the collision matrix at energy 




(5.33) 
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where the observed particle width Yr is defined from (15.251) and where the observed 
gamma width T 7ifl is given by 

r - (£) = T^Hkr (5 - 34) 

The correction factor is thus identical for the particle and gamma widths [see (15. 6p ]. 
For a resonant process, the main part of the wave function is located at short distances, 
and the external term (I5.32p can often be neglected to a good approximation. In that 
case, the capture cross section (12.571) takes the usual Breit-Wigner form 



r .A (F) „ 2J+1 T ltR (E)T R (E) 

Ji*,M ) ~ k 2 (2/i + 1)(2 j 2 + 1} {Er _ E y + {Y R {E)/2f 



As the electromagnetic interaction is weak, it is implicitly assumed that T lR (E) <C 
T R (E). This formula is of course an approximation which assumes that (i) there is 
no background or other resonances interfering and (ii) the external contribution is 
negligible. Going beyond these two approximations can be done by using the more 
general formulas (I5.29P and (15. 32p . Notice that the relative roles of and U ext depend 
on energy. When E becomes small, the internal part of the initial state becomes smaller 
and smaller, and the importance of the external contribution increases. The contribution 
of the external term also depends on the binding energy of the final state. If Ef is small, 
the asymptotic decrease of the bound-state wave function is slow and the external matrix 
element (I5.32j) may be important. 

In the single-channel approximation, the internal and external components can be 
combined, which yields a slightly different definition for the gamma width [17] . In |123] , 
the term involving the collision matrix in (I5.32j) is recast with the internal contribution. 
This provides modified electromagnetic matrix elements. As long as the external capture 
is negligible, all definitions of the gamma width are equivalent. 



5.7.3. Application to 12 C(p,j) 13 N 

As discussed in section 4.12, the 12 C(p,7) 13 N S-factor at low energies is essentially 
determined by the properties of the l/2 + (I = 0) resonance at E R = 0.42 MeV in 
13 N \10S\ [TTj . As usual in nuclear astrophysics the main issue is to extrapolate the 
available data down to stellar energies (around 24 keV at the typical temperature 
1.5 x 10 7 K). We use the phenomenological i?-matrix approach with the single-pole 
approximation. In that case, four parameters are to be considered: the energy, proton 
widtf T R , gamma widths r 7 , R {E R ) of the resonance, and the ANC C\ f of the 13 N ground 
state (J/ = l/2~, If = 1). Since we are dealing with very low scattering energies and 
since 13 N is not strongly bound (Ef = —1.94 MeV), the ANC should be included. The 
fit is performed using (I2.57p . (I5.32p . (I5.33P at a channel radius of a = 5 fm, and the 
resulting observed parameters are given in Table [TH 

Figure [TTJ presents the .R-matrix 5* factor compared with experimental data. With 
the single-pole approximation the fit is not perfect, in particular above the resonance, 
where the i?-matrix calculation slightly overestimates the data. This was already 
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Table 14. Observed i?-matrix parameters at the resonance energy En for the 
12 C(p, 7) 13 N reaction. 

E R (MeV) T R (keV) T y>R (eV) C x (fnT 1 ^ 
0.415 31 0.4 1.1 



observed in Fig. [121 with the calculable approach. This problem has been addressed 
by Barker and Ferdous [T7] who showed that an excellent fit of the data requires at 
least two poles in the i?-matrix expansion. In addition to the total S factor, we also 
present the internal and external contributions independently. This analysis is done for 
a = 5 fm, but also for a = 6 fm, where the same observed parameters are used. The 
internal part (dashed line) corresponds to the Breit-Wigner approximation (I5.35P ; it is 
almost insensitive to the choice of the channel radius. On the contrary, the external part 
(dotted lines) does depend on a. Its influence near the resonance energy is weak but, as 
expected, it increases at low energies. Neglecting the external term when extrapolating 
down to stellar energies would provide a strong underestimation of the S factor. 




0.1 0.2 0.3 0.4 0.5 0.6 

E (MeV) 



Figure 17. 12 C(p, 7) 13 N S'-factor computed with the parameters of Table [141 at a = 5 
and 6 fm. The dashed line corresponds to the internal contribution (a=5 fm and a = 6 
fm arc indistinguishable), and the dotted lines to the external term. 
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6. Recent applications of the J?-matrix method 

6.1. Introduction 

In section 4 we gave simple examples of the calculable i?-matrix method, in order to 
illustrate the theoretical framework with applications which can be easily reproduced 
by the reader. However, in most cases, alternative methods, sometimes simpler, are 
available. 

In this section, we present more ambitious applications of the i?-matrix theory in 
nuclear physics. The first deals with microscopic cluster models [96J, where the relative 
motion between the colliding particles is not given by a potential, but by a nucleon- 
nucleon interaction. The i?-matrix theory is also very efficient to solve coupled-channel 
problems [124] . In various models, the Schrodinger equation is reduced to a system of 
coupled differential equations. This can be solved, for example, with the generalized 
Numerov algorithm, but this method looses stability when the size of the system 
increases [125] . In that case, the i?-matrix theory provides an efficient alternative, in 
particular when it is associated with the Lagrange-mesh method [93]. Two applications 
concerning the three-body continuum are presented: (i) the Continuum Discretized 
Coupled Channel (CDCC) method [1261 187] . and (ii) the three-body hyperspherical 
formalism [127[ 1128} IS3"] . Recent applications of the i?-matrix theory in atomic physics 
can be found, for example, in [72} [73], 11291 1130] . 

6.2. Microscopic cluster models 

6.2.1. General presentation 

In nuclear physics, a microscopic theory is based on a description of all nucleons taking 
full account of antisymmetrization and derived from interactions between nucleons. The 
Hamiltonian (12.251) reads 

H = £S-+ E Yak-Ten, (6-1) 

i=l N i>j=l i>j>k=l 

where A is the nucleon number, p i is the momentum of nucleon i and Vy and 
Vijk are two- and three-nucleon interactions. In (16. ip . subtracting the center-of-mass 
energy T cm guarantees that the wave function is free of spurious cm. components 
[131] . The two-body interactions Vij involve a nuclear term with spin-orbit, tensor 
and other components, and the Coulomb interaction. The three-body interactions V^k 
are necessary to explain the binding energies of the 3 H, 3 He and 4 He nuclei. Realistic 
interactions are derived from field theories and partly fitted on properties of the nucleon- 
nucleon system. For small nucleon numbers (i.e. A < 4), different techniques exist to find 
numerically exact solutions of the Schrodinger equation [132] . Few-body calculations can 
be performed by using various realistic interactions and compare well with experiment. 
For heavier systems, ab initio calculations [133] 1134] become available. Different variants 
are being developed, but are currently limited to A m 12, and their application to 
continuum states is quite difficult |135[ 1136] . 
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The microscopic cluster approach is based on an assumed cluster structure, i.e. on 
the existence of correlated subsystems in the fully antisymmetric wave function of the 
A-nucleon system [96] . The microscopic cluster model provides a unified framework for 
the description of nuclear spectroscopy and of nuclear reactions. The cluster assumption 
allows the application of a microscopic theory to heavier systems (typically up to 
A ~ 20 — 24), but requires the use of effective nucleon-nucleon two-body interactions 
(see, for example, [1371 I138j ) which are adapted to the cluster approximation. We 
use here the Minnesota effective interaction |138] which does not include tensor forces, 
but simulates their contribution in the binding energy of the deuteron by the central 
term. Three-body effects are in general also approximately simulated in these effective 
interactions. 



6.2.2. The Resonating Group Method 

The Resonating Group Method (RGM) assumes that the wave function can be 
expressed in terms of normed cluster wave functions <pi involving Aj, among the 
A nucleons and depending on translation-invariant internal coordinates, and of an 
unknown wave function gi(r) for their relative motion. Using the isospin formalism, 
the A-nucleon approximate wave function reads 

% m = jf^AfafogiirWrm, (6.2) 

where r = (r, Q) is the relative coordinate between the cm. of the clusters, and A is 
the A-nucleon antisymmetrization projector 

^=^E(-im, (6-3) 
A - P =i 

where the operator P p performs the permutation p among A particles. For simplicity, 
we assume a two-cluster structure, with zero-spin clusters [see (12.321) for a multichannel 
definition]. More general presentations can be found in [961 EH1 EH I139[ 1140] . The 
antisymmetric internal wave functions <pi and <p2 ar e normalized to unity and defined in 
the harmonic-oscillator shell model, with a common oscillator parameter b [96]. First 
applications were limited to s-shell clusters [35] , but progressively, heavier systems were 
considered [T4"T] [T42"] . 

The unknown relative function gi(r) is obtained by solving the Schrodinger equation 

HV lm = EV lm , (6.4) 

with Hamiltonian (16.11) and approximate wave function (16. 2p . By projecting the 
Schrodinger equation over (0i</>2YJ m (f2)|, this technique provides an integro-differential 
equation involving local and non-local potentials [139J 

(T r + V D (r) - E) 9l {r) + J™ (jC Hl (r,r>) - EK Nl (r,rj)g l (r')dr' = 0, (6.5) 

where T r is the relative kinetic energy operator ( 12. 39ft . Vd{t) is the direct potential 
(I2.40p . and /C^(r, r') and K.Ni(r,r') are the Hamiltonian and overlap exchange kernels, 
respectively. 
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At large distances, the Hamiltonian (16. ip can be written as 

H — > Hi + H 2 + T r + V c (r), (6.6) 

1 — >oo 

where Hi and H 2 are the internal Hamiltonians of the clusters, similar to (16. ip . and 
Vc(t) is the Coulomb potential between charges Z\e and Z 2 e. The internal energies E{ 
of clusters 2 = 1,2 are defined by the variational expressions 

E % = (falHilfa). (6.7) 

In parallel, the wave function (16. 2p tends to 

^ ^i&fl WW (6-8) 

since the antisymmetrization operator ^4 acts at short distances only. For large r values, 
the radial wave function rgi(r) is given by the Coulomb equation (j2.4p . Although it is 
also suitable for bound states, the ansatz (I6.2p of the total wave function is therefore 
well adapted to the treatment of scattering states. 

The main problem of the RGM is not to solve (16. 5p . This can be done by standard 
techniques, or by the R- matrix method (12] . Recent calculations, using realistic nucleon- 
nucleon interactions, have been performed on the 3 He+p and 3 H+n scattering |143] and 
on the a+nucleon scattering |136j . In these references, the RGM equation (16. 5p is 
extended to a multichannel generalization. In |136j . it is solved on a Lagrange mesh. 
The drawback of the RGM is that the calculation of the overlap and Hamiltonian kernels 
is not systematic. The reason is that the relative coordinate r and the cluster internal 
coordinates in (16. 2p are modified in different ways by the terms of the antisymmetrizer 
A. This method requires heavy analytical calculations (96J I139j (see also (321 1144] for 
examples of kernels). This problem is simplified by using the Generator Coordinate 
Method (GCM), described in the next subsection. 

6.2.3. The Generator Coordinate Method 

In the GCM, the relative wave function gi{r) is expanded over projected Gaussian 
functions |145] as 



9i(r) = J fi(R)Ti(r,R)dR, (6.9) 
where R is the generator coordinate and T; (r, R) is defined as 
r( m f/M 3/4 ( r 2 + R 2 \. ffx rR\ 

r,(r ' fl) = wJ exp ^°^^JH^J- (6 - 10) 

In this equation, /xq is the reduced mass in units of the nucleon mass, and ii(x) = 



7r/2x// + i/ 2 (x), Ii+i/2(x) being the modified spherical Bessel function of the first kind 
9J. The calculation of gi(r) is therefore replaced by the calculation of the generator 
function fi(R). 

Inserting (I6.9p in the RGM definition (16. 2p provides 

tf, m = / fi(R)$ lm (R)dR, (6.11) 
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where 

*UR) = j^A^r^R^rin). (6.12) 

After multiplication by an appropriate factor depending on the cm. coordinate of the A 
nucleons, the basis function $; m (i?) can be expressed as a projected Slater determinant 
provided that the oscillator parameters of the clusters are identical |145] . The Slater 
determinant is defined from A\ and A2 shell-model orbitals centered at — A2R/A and 
A\Rj A for the first and second cluster, respectively. This property is well adapted to 
systematic and numerical calculations since it involves matrix elements of single-particle 
orbitals only. Well-known techniques exist to determine matrix elements between Slater 
determinants from single-particle matrix elements [146^ 1147] . They allow a rather simple 
extension of the cluster model to heavy systems |148] and to multichannel problems |149] . 
The projection over angular momentum / can be performed numerically |146j . 

In practice, the integral in (16. 9p is replaced by a finite sum over a set of values 
R n of the generator coordinate. This means that, at large distances r, the radial 
wave function gi(r) presents a Gaussian behaviour, not consistent with the physical 
asymptotic behaviour. This problem can be addressed by using the microscopic R- 
matrix method [SB]. The wave function is approximated in the internal region by a 
discretized version of (16.111) as 

N 

*ta = J2fl( R n)®lm(Rn). (6.13) 
n=l 

In the external region, it is approximated by the asymptotic expression ( \6.8\\ as 

Kt = 0i«»*T(^) (6.14) 

where the external radial function rg ; cxt (r) is a linear combination of Coulomb functions, 
as in (13.31) . 

The application of the i?-matrix method to the GCM is straightforward and follows 
the method detailed in section 13.21 The generalized matrix C defined in (13.141) involves 
matrix elements of the Hamiltonian in the internal region only. This is achieved by 
subtracting the external contributions [SSj- By definition of the channel radius a, 
antisymmetrization effects and the nuclear interaction are negligible in the external 
region. The relevant matrix elements are therefore given by 

($l{Rn)\MRn>))i a t = (MRn)\MRn>)) ~ J ^(r, R^T^T, R n >)r 2 dr, 

($Ki?„)|#|$^n')>int = ($ z ( J R n )|i/|$,( J R n ,) a ) 

roo 

- T l (r,R n )(T r + V c (r) + E 1 + E 2 )T l (r,R nl )r 2 dr, (6.15) 

J a 

where the first terms in the r.h.s. are matrix elements over the whole space, involving 
Slater determinants. The second terms represent the external contributions of the basis 
functions (I6.10p and can easily be computed numerically. Then the i?-matrix and the 
associated collision matrix are obtained as in section 3.2. Similarly, the collision matrix 
should not depend on the choice of the channel radius a, provided it is large enough 
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to make the nuclear interaction and the antisymmetrization effects negligible in the 
external region. A generalization to multichannel systems can be found in [36| [37]. 

6.2.4- Applications: a + a and 12 C+p 

We present three typical applications of the GCM associated with the microscopic R- 
matrix method. The first deals with the well known a + a phase shifts (see section |4~4|) . 
Then we compute the 12 C+p elastic cross section, as well as the 12 C(p,7) 13 N S factor 
at low energies. In all cases we use the Minnesota (MN) effective interaction |138j as 
central nucleon-nucleon force. The Minnesota potential provides the correct binding 
energy of the deuteron (without tensor force) and reproduces fairly well some properties 
of nucleon-nucleon scattering. It involves the admixture parameter u whose standard 
value is u — 1, but which can be slightly modified to fit important physical quantities, 
such as the energy of a bound state or of a resonance. For the 12 C+p system, a zero- 
range spin-orbit force (with amplitude S ) is added |150j . 

The a + a system is described by cluster wave functions (pi defined from four Os 
harmonic-oscillator orbitals with an oscillator parameter b = 1.36 fm for the two protons 
and the two neutrons. With the MN force, the binding energy of the a particle is 24.28 
MeV, which is smaller than the experimental value 28.30 MeV. This difference does 
not play an important role, since all theoretical energies are defined with respect to the 
a + a threshold. The admixture parameter u is taken as u = 0.94687, as recommended 
in |144j . This value provides an excellent description of the a + a phase shifts in a wide 
energy range. We use N = 10 basis functions in ( I6.13P with R n values ranging from 0.8 
fm to 8 fm by steps of 0.8 fm. In Table [T5"l we give the + , 2 + , 4 + phase shifts at typical 
energies, and for various conditions of the calculation. The channel radius a is taken 
as a = 6.4 fm or a = 7.2 fm, and N = 9 or 10 are considered. In all cases, the phase 
shifts are very stable when the conditions are changed. They can be obtained with an 
accuracy better than 0.1°. Fig. [TBI shows the phase shifts as a function of energy. It is 
known that a cluster model is well adapted to the a + a system since the a particle has 
a large binding energy and since the first open threshold ( 7 Li+p) is near 17 MeV. The 
GCM phase shifts are therefore in very good agreement with experiment. The same 
results can be obtained from the RGM equation ( 16. 5ft [144] . 

The second application deals with the 12 C+p elastic scattering. The 12 C internal 
wave function (b = 1.65 fm) is described in thep shell limited to zero angular momentum 
and isospin. This provides four independent shell-model states defined from a linear 
combination of Slater determinants [152] . The calculation is therefore performed with 
four channels, obtained from the diagonalization of the 12 C basis. The first eigenstate 
corresponds to the ground state, and the three additional eigenvalues are considered as 
virtual excitations, which simulate the distortion of 12 C in the 12 C+p system. 

The spin-orbit amplitude is fixed as So = 36.3 MeV.fm 5 , which allows to reproduce 
the l/2~ and 3/2~ binding energies of 13 N with u = 0.77. This value is used for negative- 
parity partial waves. In positive parity, u = 0.998 reproduces the experimental energy 
of the l/2 + resonance. The GCM widths of the l/2 + and 3/2~ resonances are 43 keV 
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Table 15. Microscopic a + a phase shifts (in degrees) for different conditions of 
calculation. 
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Figure 18. a + a phase shifts with the microscopic cluster model. Experimental data 
are taken from |151j . 



and 99 keV, somewhat larger than the experimental values (31.7 ±0.8 keV and 62 ± 4 
keV, respectively [9~4"]). 

The 12 C+p excitation functions at 9 = 89.1° and 9 = 146.9° are shown in Fig. [191 
Below 1.5 MeV, the agreement with experiment is quite good. However the resonant 
structure near 1.6 MeV is not well described since the experimental data involve a 5/2 + 
resonance, which cannot be described with a simple cluster structure. This problem 
often occurs in microscopic cluster models, where some experimental states do not have 
a cluster structure, and are therefore not present in the model (see, for example, |153] 
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Figure 19. 12 C+p excitation functions at two cm. angles. The data are from |120j . 



for the 16 0+p system). Of course the phenomenological approach (see Fig. H~5l) provides 
a better description, but the parameters are adjusted to the data, and all resonances 
are included from the very beginning. 

Our third application of the GCM is the 12 C(p, 7) 13 N cross section at astrophysical 
energies. The matrix elements of the El operator are determined as explained in 
section 15.71 (see [37] for detail). The capture cross section requires the 13 N ground 
state wave function, as well as 12 C+p scattering states. Only s waves, corresponding to 
the l/2 + resonant partial wave are included in the scattering state. Calculations with d 
waves show that this component is of the order of 15% at 1 MeV, but becomes negligible 
as soon as the energy decreases, because of the higher centrifugal barrier. 

The resulting S factor (12.551) is shown in Fig. [20j As expected from the 
overestimation of the proton width, the peak near the l/2 + resonance at 0.42 MeV is 
slightly too broad. This type of calculation cannot be expected to perfectly reproduce 
the data. However the S factor is quite satisfactory considering the fact that no 
parameter is fitted to capture data. The model has thus a predictive power which 
is very useful for capture reactions that have not been measured yet or can not be 
measured. The GCM has been used to study many reactions of astrophysical interest, 
to compute either the cross section or the properties of low-energy resonances (see for 
example [154] ). 
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Figure 20. 12 C(p, 7) 13 N GCM S-factor. The experimental data are from [107] (open 
circles) and |108j (closed circles). 



6.3. The Continuum Discretized Coupled Channel (CDCC) method 
6.3.1. Formalism 

The purpose of the CDCC method is to determine, as accurately as possible, the 
scattering and dissociation cross sections of a nucleus which can be easily broken up 
in the nuclear or Coulomb field of a target. The final states may thus involve three 
particles: the target and the fragments of the projectile. The relative motion of these 
fragments is described by approximate continuum wave functions at discrete energies. 
The CDCC method was suggested by Rawitscher [126] and first applied to deuteron + 
nucleus elastic scattering and breakup reactions. It was then extensively developed and 
used by several groups [871 1125] . Its interest has been still revived by the availability of 
radioactive beams of weakly bound nuclei dissociating into three fragments, such as 6 He 
whose first dissociation channel is a+n+n [1551 1156] , Although variants of the CDCC 
method also exist in atomic physics [157] . we focus here on applications in nuclear physics 
with two-body projectiles. To simplify the presentation, we assume a spin zero for the 
constituents of the projectile, and for the target t. As usual in CDCC calculations, the 
internal structure of the three particles is neglected. 

Let us consider the coordinate system of Fig. [2TJ R is the internal coordinate of the 
projectile and r is the coordinate of the relative motion between projectile and target. 
The three-body Hamiltonian is given by 

H = H + T r + V tl (r + ^R) + V t2 (r - ^-R) , (6.16) 

where H Q is the two-body Hamiltonian of the projectile 

H = T R + V 12 (R). (6.17) 

In general, potential Vu(R) associated with the projectile is real, whereas the 
interactions Vu between the fragments i and the target t are derived from the optical 
model, and thus complex. In a schematic notation, the wave function associated with 
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(16 . 161) can be expanded as 

tt(r, R)=Y, ®B(R)u B (r) + / <5> k (R)u k (r)dk, (6.18) 

where B denotes the bound states of the projectile, and 4>k(R) are two-body scattering 
wave functions with wave number k. The first term represents the elastic and inelastic 
channels, and the second term is associated with the breakup contribution. 

In practice, two methods are available to perform the continuum discretization, i.e. 
discretize the integral over k. In the "pseudo-state" approach, it is replaced by a sum 
over square-integrable positive-energy eigenstates of Hamiltonian (16.171) . The projectile 
Hamiltonian H is diagonalized over a finite basis, yielding the square-integrable radial 
functions $f (R) at energies E^, 

H <S>f(R)Y L M ({l R ) = E^f(R)Y L M (Q R ). (6.19) 

These functions are associated with bound states (i = B,Ei < 0), or represent square- 
integrable approximations of continuum wave functions (Ei > 0). 

The alternative is to replace the integral over k by averages of exact scattering 
states over a range of energies ( "bin" method) [87] . This approach also provides square- 
integrable basis functions. As far as the applicability of the i?-matrix is concerned, both 
methods are treated in the same way. We use here the pseudo-state method. 

1 

— — o, 

2 

Figure 21. Coordinate system used for CDCC three-body calculations. 




The total wave function (I6.18P is then rewritten, for an angular momentum J and 
parity n = (— l) l+L , as 

* JM «(r,R) =E Y ii M (nr,n R )<S>t(R)u£(r), (6.20) 

ILi 

where J results from the coupling of orbital momenta I and L, and 



JM 



The relative wave functions uf^ir) are given by a set of coupled equations 



h 2 d 2 

2/x \dr 2 



1(1 



+ E 



E 



(6.21) 



(6.22) 



where the channel index c stands for (ILi). Of course, the sum over L must be truncated 
at some value L max . The sum over the pseudo-states i is limited by the number of basis 
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states and can be reduced further by eliminating states above a maximum energy E max . 
The CDCC problem is therefore equivalent to a system of coupled equations where the 
potentials V^f(r) are given by 



V, 



Jn 

<:<:' 



<X^(nr,n R )9f(R)\V a (r + -±R) + V t2 (r - ^r) \Y™(n r , n R )Qjf(R)). (6.23) 

This matrix element represents a 5-dimensional integral over (Q r , Q R , R). In practice 
the potentials are expanded into multipoles as 

vJr + ^R) + vjr - ^R) = £F A (r, R)P X (cos 9 Rr ), (6.24) 



A J \ A 

where 9 Rr is the angle between R and r, and P\[x) a Legendre polynomial. In practice 
the number of A values is limited by angular-momentum couplings. The four angular 
integrals in (I6.23P are performed analytically, whereas the integration over R may require 
a numerical approach. 

The system (16.221) can be solved by various methods |124] , in particular with the R- 
matrix formalism. We expand the radial wave functions u J c ™{f) over Lagrange functions 
( 14. 7j) . and the calculation of potential (16.231) is therefore limited to the mesh points as in 
(I4.12p . Further simplifications are possible by describing the projectile wave functions 
$f (R) [see (16.191) ] in a Lagrange basis as well. In this well known two-body problem 
|93j . the wave function is expanded in a basis differing from ( 14 .7p by the fact that it 
is constructed from Laguerre polynomials. The matrix elements involving $f (R) [as, 
for example, in the potentials (16.231) ] are determined in a fast and accurate way. The 
extension to three-body projectiles |155l 1156] can also be considered in this approach. 
The calculations are much more time-consuming since the projectile wave functions 
depend on two radial coordinates, and are more difficult to handle. Consequently, 
even though the formulation is similar, the calculation of the potential matrix elements 
(I6.23P raises important numerical difficulties, which could be addressed by the present 
technique. 



6.3.2. Application to the d+ 58 Ni elastic scattering 

The CDCC theory, associated with the i?-matrix method, is applied to the elastic 
scattering of deuterons on 58 Ni at E\ a ^ = 80 MeV. This collision has been widely covered 
in the literature [87 ] 1158] . The p + n and nucleon- 58 Ni interactions are chosen as in these 
references. 

The first step is to determine the deuteron ground state, and the p + n pseudo- 
states, from (16.191) . These wave functions are expanded over a Lagrange-Laguerre basis, 
involving a scaling parameter h which is adapted to the size of the system [93] . Typical 
values are h ~ 0.3 — 0.4 fm with 15 ~ 20 basis functions. We include partial waves 
L = 0,2,4. In a second step, the potentials V cc /(r) [see (16.221) ] are computed at the 
mesh points of the d+ 58 Ni relative motion. These mesh points are defined from the 
channel radius a and from the size of the basis N [see (I4.8p ], 



The R-matrix Theory 



77 



Table 16. Amplitude ?y /7r and phase shift <5 j7r for elastic d+ 58 Ni elastic scattering 
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In Table [16] and Fig. [221 we present the elastic part of the collision matrix 

jjJn = ^ eX p(2 Z( ^) 

for J 71 " = + and 17". Values for J 71 = 17" can be compared with the literature |158j . 
The amplitude r\ and phase shift 8 are computed for different channel radii and numbers 
of wave functions. For the pseudo-states, we truncate at L max = 4 and E msix = 50 
MeV. The number of pseudo-states depends on the size of the basis (here, typically 
their number is about 20). We have checked that changing the size of the basis does 
not affect the phase shifts. Table [TBI shows that, for high partial waves (J* = 17"), the 
stability with the numerical conditions is virtually perfect. The amplitude is in excellent 
agreement with |158j . but the phase shift is slightly different (0.5°). For low partial waves 
(J n = + ), the stability is still acceptable but higher N values are necessary. This can 
be understood by the fact that the internal part is more and more important as J 
decreases. 

0.66 
0.64 
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V 0.62 
0.60 



0.58 

20 40 60 

E max (MeV) 

Figure 22. Amplitude i] 17 ~ for different L max values (labels) as a function of the 
truncation energy E max . 

Fig. [22] displays the amplitude t] 17 as a function of the p + n truncation energy 
E max , and for different L max values. The convergence with respect to E m ^ is reached 
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near E max » 40 MeV, which corresponds to /c max ~ 1 fm _1 . This result agrees with the 
conclusion of |158j . From the figure, it is clear that L = pseudo-states are not sufficient 
to provide accurate values. However, a very good convergence is already obtained with 

^max 2. 

The elastic cross section is presented in Fig. [231 and compared with experimental 
data (quoted in |159j ). A first calculation, referred to as "no-breakup" approximation, 
is done by including only the ground state of the deuteron (dotted line). This 
approximation clearly overestimates the data above 30°. With L = deuteron pseudo- 
states, the agreement in the region 30° — 50° is significantly improved. The present 
result is very close to the CDCC calculation of |159j . The curves with L max = 2 and 
L max = 4 are indistinguishable at the scale of the figure. 




Figure 23. d+ 58 Ni elastic cross section relative to the Rutherford cross section at 
E lah = 80 McV. Solid lines correspond to E max = 40 MeV and to L max = 0, 4. The 
dotted line corresponds to the no-breakup approximation. Experimental data are taken 
from [155]. 



From this example, it turns out that the i?-matrix method is a very efficient tool 
to solve the CDCC equations. By using Lagrange functions, the calculation of the 
coupling potentials V cc /(r) is very fast, since only values at the mesh points are required 
(typically N ~ 30 — 40). Consequently, the main part of the computer time is devoted 
to the inversion of the (complex) matrix C [see fl3.15p ]. Fast techniques are available for 
matrix inversion, and the total computer time is therefore very short (typically a few 
seconds for the angular distribution of Fig. 1231) . This opens encouraging perspectives 
for CDCC calculations involving three-body projectiles which require large computer 
times with standard techniques [1551 1156] . 



The R-matrix Theory 



79 



6-4- Three-body continuum states 

6.4-1- Hyperspherical formalism 

The separation energy being very low in exotic nuclei, a precise treatment of 
the continuum is necessary for the description of reactions leading to a three-body 
dissociation. The i?-matrix method provides an efficient way to treat three-body 
continuum states |128[ 183] . We use the hyperspherical formalism, well adapted to three- 
body systems [1271 ESI I160j . This model is well known, and is just briefly outlined 
here. 

In a three-body system, the Hamiltonian is defined, after removal of the cm. 
motion, as 

3 3 

H = ^T,-T cm + E V u( r j ■'•>)■ (6.25) 

i=l «>i=i 

where 77 are the space coordinates of the particles, Tj their kinetic energy and Vij the 
interaction between particles i and j. The scaled Jacobi coordinates [x = (x,Q x ) and 
y = (y, fly)} are defined from the coordinates r i: and provide the hyperradius p and 
hyperangle a as 

P = \j% 2 + y 2 , 

a = arctan(y/x). (6.26) 

In this coordinate system, the 3-body kinetic energy involves the operator K 2 which 
generalizes the concept of angular momentum in 2-body systems. It commutes with 
l 2 x and I 2 and their common eigenfunctions (^5) are known analytically (see 
[88] for details). The eigenvalue of K 2 is K[K + 4) where the integer K is the 
hypermomentum quantum number. In these expressions, f2 5 = (Q x , Q y , a) and 7 stands 
for 7 = (l x , l y , L, S) where (/ x , l y ) are the angular momenta associated with (x, y), L and 
5* are the total orbital momentum and spin, respectively. The total angular momentum 
J results from the coupling of L and S, and the parity is given by it = (— 1) . 

The wave function of partial wave Jir associated with Hamiltonian (I6.25P is then 
expanded over hyperspherical harmonics as 

* JMn (pM = P" 5/2 EX^(P) y^m, (6.27) 

-yK 

where the radial functions X^k{p) have to be determined. The Schrodinger equation is 
then reduced to a system of coupled differential equations 

h 2 f<P_ _ (K + 3/2) (if + 5/2) 
2m N \dp 2 p 2 

+ E ^w>) *~Mp) = 0, (6-28) 
where the potentials matrix elements are defined as 

^Wp) = OtfZWl E ViArj - ri)|3^(n 5 )>. (6-29) 

i>j=i 
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The integral over Vt x and Vt y are performed analytically, whereas a numerical quadrature 
is used for the integral over the hyperangle a. With the Raynal-Revai coefficients 
|161l 1162] the evaluation of (16.291) is rather easy. A truncation must be done in the 
summation over K\ this provides a maximum K value, denoted as K max . The number 
of components in ( I6.27P increases rapidly when K max increases. 

This formalism has been extensively applied to three-body bound states [88]. In 
|163j . we have addressed this problem by expanding the radial wave functions over 
Lagrange functions. This approach provides a very fast method to evaluate matrix 
elements flBT29l . 

One of the main issues associated with nuclear three-body problems is the presence 
of forbidden states in the nucleus- nucleus interaction [961 EH]- These two-body forbidden 
states introduce spurious states in the three-body problem and should be removed. 
This problem has been discussed in detail by Thompson et al. |128] . and is usually 
solved, either by introducing a projector over the forbidden states |164j . or by using 
supersymmetric potentials |165j . 



6.4-2. The R-matrix method for three-body states 

The treatment of three-body slates, with exact three-body asymptotic conditions, is 
more recent |128l 183] . In the i?-matrix method, the solutions of the system (16.281) are 
split in two regions, 

N 

xS,int(p) = E^^(p), (6-30) 
1=1 

for p < a where (fi(p) are basis functions and 

X^,ext(p) = C£ [H- K (kp)5 w 5 KK , - U^, K ,H+ K (kp)] (6.31) 

for p > a. In this equation, is a normalization coefficient, U JlT is the three-body 
collision matrix and the incoming and outgoing functions H^ K {x) are defined as 

H jk( x ) = ±* ( ~2 J I j k+2(x) ± iY K+2 (x)] , 

where J n (x) and Y n (x) are Bessel functions of first and second kind, respectively. We 
assume here systems without two-body Coulomb interaction. 

As in previous applications, we choose Lagrange functions for the basis states (fii(p) 
[83J. The i?-matrix method is then used to solve the coupled system (16. 28ft . Formally 
this is equivalent to the CDCC system [see (I6.22p ]. although the external wave function 
involves Bessel functions. Another difference arises from the long range of the three-body 
potential (I6.29p . As shown in [83J, this potential behaves as 

^WP) (6.32) 

even with short-range two-body interactions. This arises from the definition of the 
hyperspherical coordinates. Even for large p values, two particles can still be close 
to each other and interact strongly. Constants V^/y ^ can take rather large values 
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(examples are given in [83]). For this reason, the R- matrix radius should take large 
values (typically a ~ 200 — 300 fm) to ensure that potential ( 16 .32 j) is negligible compared 
with the centrifugal term in ( 16.281) . In those conditions, the propagation techniques 
described in section 13.91 are necessary to avoid huge basis sizes. 

6.4-3. Application to a + n + n three-body scattering 

The 6 He nucleus is an ideal test case for three-body continuum states and has been 
considered in previous studies |128[ [83] . Accurate a + n and n + n interactions exist in 
the literature |166l I138j . The a + n subsystem presents one forbidden slate for / = 0, 
which is removed by using a pair of supersymmetric transformations |165] . Details can 
be found in [S3]- 

Figure EU illustrates the need for propagation techniques. Dotted lines correspond 
to channel radii a = 20 fm and a = 30 fm, and are compared with the exact values 
(solid line, a = 250 fm) obtained by propagation. As expected from the long range of 
the three-body potential (16.291) . large values for the channel radius are necessary. As 
soon as this condition is satisfied, the .R-matrix phase shifts are very stable against 
variations of a. 

In Fig. [251 we display the + eigenphases as a function of K ms , x . The eigenphases 
are obtained after diagonalization of the collision matrix [see ( 12. 37ft ] . In each case 
we select the eigenphase with the dominant resonant structure. Above 4 MeV, a 
fair convergence is obtained, but high hypermomenta are necessary near 2 MeV. The 
+ , 1~, 2 + eigenphases are displayed in Fig. |26]with i^ max = 24, 19, and 16, respectively. 
The 2 + phase shift presents an experimentally well known narrow resonance at low 
energies. For the + and 1~ partial waves, the calculation shows a broad structure near 
1.5 MeV. The existence of three-body resonances at low energies, and in particular for 
J 71 = 1~, is still an open debate (see, for example, the discussion in |167j ). from the 
experimental as well as from the theoretical viewpoints. 

In the future, this formalism could be applied to the 3a system. Current 
experimental |168] and theoretical [1691 1170} 1171] results are rather controversial 
concerning the existence (or non-existence) of broad + and 2 + resonances above the 
3a threshold. This issue is crucial in nuclear astrophysics, since the Hoyle state (0^) 
in 12 C is known to be the main resonance in Helium burning. If the properties of this 
resonance are well known, the 12 C level scheme above this resonance is still not clear. 
In this context, the calculation of 3a phase shifts would help clarifying the situation. 
However, existing a + a local potentials do not provide a satisfactory description of the 
12 C spectroscopy [172j . Non-local potentials are more promising [17211173] . but although 
bound states can be easily investigated with Lagrange meshes [83], their application to 
three-body continuum states remains a challenge for the future. 
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Figure 24. a+n+n phase shifts ( J v = + ) for channel radii a = 20 fm (iV = 20) and 
a = 30 fm (AT = 30) without propagation (dashed lines). Solid lines are obtained with 
propagation up to a = 250 fm (from [SSj)- 





Figure 25. Energy dependence of a+n+n cigenphascs (J 71 ' = + ) for different K, 
values (from |83j). 



The R-matrix Theory 



S3 



300 
240 
180 

D) 
CD 

to 120 
60 




6 He 








--"-""^ 2 + 




1 - - 










J 







4 6 
E (MeV) 



10 



Figure 26. Eigenphases of 6 Hc for different J values (from |83j). 
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7. Conclusion 

The R matrix method was born sixty years ago with the important but rather limited 
goal of describing resonances in nuclear reactions. Today it has evolved into powerful 
tools covering several subfields of atomic, molecular and nuclear collisions. The literature 
is so enormous that it is not possible to master it and to cover it in a single review. 

We have emphasized a fact that is often unknown to R matrix practitioners: two 
variants exist. The phenomenological R matrix remains close to the original spirit 
and is very much used in nuclear physics to parametrize low-energy cross sections. Its 
main merits are that all parameters are real and that they have a physical meaning. 
Although resonances often play a crucial role in these parametrizations, non-resonant 
cross sections are accurately described as well. The calculable R matrix is an efficient 
technique to solve the stationary or time-dependent Schrodinger equation in various 
situations as well as its relativistic extensions. It underwent most of its developments 
in atomic physics but we have shown in section [6] that it can also be useful in nuclear 
physics. 

Because of the variety and complexity of its applications, the i?-matrix theory has 
been and is still sometimes misunderstood or misjudged for different reasons. A first 
often-made criticism concerns the role of the channel radius. The i?-matrix formulas 
depend on this radius which has no strict physical meaning. This criticism can not be 
addressed to the calculable R matrix where the independence of the physical results 
on this radius, provided that it is large enough, is a useful validity test. In the 
phenomenological .R-matrix however, the channel radius is indeed a parameter whose 
value is disputable. This arises from the truncation of the R matrix to a small number 
of poles. In spite of this truncation, this approximation often provides excellent fits 
to the data. The standard option is to optimize the channel radius together with the 
other parameters. This radius should however be larger than the sum of the radii of 
the colliding nuclei. Its value should always be mentioned because the other parameters 
depend on its choice. 

Another criticism deals with a reputation of poor convergence of the calculable R 
matrix. As we have shown, this reputation is undeserved. Its origin lies in the choice of 
a common logarithmic derivative for the basis states in the founding papers. While this 
choice is acceptable (although with some discontinuity drawback) for an infinite basis, 
it leads to inaccuracies when the basis is truncated. The introduction of the Bloch 
operator has opened a way to the use of finite bases providing a sufficient variety of 
behaviours at the boundary for which this problem disappears, a fact not yet known 
enough. Modern i?-matrix codes can employ different types of such bases. They reach 
an excellent convergence and do not require the use of the Buttle correction. 

The calculable R matrix provides bound-state and scattering-state wave functions 
that can be used in a variety of applications with sometimes technical complications due 
to the existence of two regions to define the wave function. In atomic, molecular and 
nuclear physics, the challenge is now to reach the same level of accuracy for processes 



The R-matrix Theory 



85 



with more than two unbound particles in the final states. Advances have been made 
for double ionization by a photon in atomic and molecular physics and for the breakup 
of two- and three-body halo nuclei in nuclear physics. Progresses of i?-matrix theory in 
these directions should still be expected. 
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A. Appendix: Collision matrix and K matrix 

System (12.381) possesses N linearly independent solutions that vanish at the origin. A 
matrix solution u is obtained by putting those independent solutions as columns of a 
square matrix. By multiplication on the right by any invertible matrix, one obtains 
another matrix solution of system (I2.38H which is physically equivalent. The most 
general asymptotic expression of such a matrix generalizes (12.34)) as 

u — ► v-V 2 (I - OU)C, (Al) 

where U is the collision matrix, I and O are complex conjugate diagonal matrices 
involving incoming and outgoing Coulomb functions (I2.8P on the diagonal, and v is a 
diagonal matrix of velocities. Complex matrix C is arbitrary non singular. In ( 12.34)) . 
we have chosen C diagonal for simplicity. 

With C = i(l + U)~ 1 C where matrix C' is also non singular, an equivalent 
asymptotic form which is often used is obtained with (12. 8 j) as 

u — > v~ 1/2 (F + GK)C (A2) 

where F and G are real diagonal matrices involving regular and irregular Coulomb 
functions Fi and G\ on the diagonal. This asymptotic form is real if C' is real. It is 
then the most general asymptotic form of real solutions. 
Matrices U and K are related by 

U= (1 -iK)-\l + iK) (A3) 

or 

K = i(l-U){l + U)- 1 . (A4) 
Matrix K is real and symmetric if U is unitary and symmetric. 

B. Appendix: Proof of relation ( 13.271) 

Let B be an invertible N x N matrix and u and v be vectors with components. The 
inverse of the square matrix 

A = B + uv T (Bl) 

is given by 

A- = B-'- B "'"f B 1 "', (B2) 

1 + v T B^u V ; 



where the denominator is a scalar. A corollary of ( IB2p reads 

A- l u= . (B3) 

1 + v 1 B u 

Another corollary is the relation 

(i^A^u)- 1 = 1 + ^B^u)- 1 (B4) 
from which (I3.27P follows. 
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C. Appendix: Matrix elements for various basis functions 

Here we present the matrix elements used for different basis functions in Sect. 4. Unless 
specified otherwise the matrix elements of the kinetic energy are given for I = and in 
units of h 2 /2/i. 

(i) Sine functions 

The overlap matrix elements between basis functions (14. 3[) are given by 

faibi) = \ 8 ij- ( C1 ) 
The matrix elements for the kinetic energy are simple, 

7T 2 /. l^ 2 



((PilTofa) = — f % - - j 5ij. (C2) 

Because of property (I4.4p at the boundary, matrix elements of the Bloch operator 
£(0) vanish. The matrix elements of 1/r 2 are also analytical but involve Sine 
Integral functions. For the potential, a numerical treatment is necessary. 

(ii) Gaussian functions 
Let us define 

pa 

IkW) — I r k exp(—ur 2 ) dr 
Jo 

= 7 ((A; + l)/2,z/a 2 )/2z/( fc+1 )/ 2 , (C3) 

where 7 is the incomplete gamma function and a is implied. The overlap matrix 
elements between basis functions (14. 5 p are given by 

(ipi\ipj) = hi+2{vi + Vj), (C4) 

with Vi = 1/b 2 . For the kinetic energy, we have 

(ifi\Ti + C(0)\<pj) = AuiUjlx+tfa + Vj) 

-2(l + l)(v i + v j )hi +2 (v i + v j ) 

+ {l + l){2l + l)I 2l {v i + v j ). (C5) 

For a Gaussian potential and for the Coulomb potential, the matrix elements read 

(<Pi\ exp(-(r/r ) 2 )|v3j) = hi+2^i + + 1/rj}), (C6) 
((Pi\l/r\(pj) = hi+ii^i + Vj), (C7) 

and therefore do not require any numerical integration. Of course, other potentials 
can be considered, but the matrix elements must be, in general, obtained from a 
numerical integration. 

(iii) Lagrange functions 

Let us start with matrix elements in interval (0, a). The regularization coefficient 
n in (j4.7p is taken as n — 1. This ensures that the wave function vanishes at 
the origin and allows that the Coulomb potential is treated accurately at the 
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Gauss approximation. Using wave functions (14.91) and the corresponding Gauss 
approximation the matrix elements take a very simple form 

{<Pi\<Pi)=6v (C8) 

(<p i \V\<p j ) = V(as i )5 ij (C9) 
i.e. they only require the evaluation of the potential at the mesh points. For the 
kinetic energy, a simple calculation [36] provides for % — j 

(47V 2 + AN + 3)xi(l - n) - 6x, + 1 



and for % ^ j, 



(<Pi\To + £(0)\<Pi) 



(<^|T + £(0)|^) 



3a 2 x 2 (1 — Xi) 2 



(CIO) 



x 



N 2 + N + 1 



a 2 [0^(1 - Xi)(l - Xj)] 1 / 2 

X 2 I «3j 2 37 j 37 j 1 



l 



[Xj, Xj 



1 37 ? 1 X 7 



(Cll) 



Thanks to the Bloch operator, this matrix element is symmetric. For 1^0, the 
centrifugal term is included in the potential (1C9|) . 

Next we consider (ai, 02) intervals which are used in the propagation method. The 
basis functions ( 14. 7ft with n = are extended to 



^1/2 Pjy((2r-a 1 -a 2 )/Aa) 
r — XjAa — ai 



(C12) 
(C13) 



with Aa = 02 — 01, and the Lagrange condition becomes 

<Pi(a,i + XjAa) = (AjAa)" 1 / 2 ^-. 
The matrix elements of the potential read 

(<Pi\V\(pj) = V(ai + XiAa) (C14) 

and are still given by a simple evaluation of the potential at the mesh points. The 
matrix elements of the kinetic energy are given at the Gauss approximation by 



for i — j and 



{<Pi\T \<Pi) 



{<Pi\T \<Pj) 



3Ao 2 Xi(l 



N + N + 6 



Xi(l Xi) 
Xj(l — Xj) 2 x 2 x j I 3 x ^ x j 4 x ^ 



Aa 2 \xi(l-Xi) Xi(l - Xi)(xj - x^ 2 
for i 7^ j. The matrix elements of the Bloch operators read 



and 



(Vil^oalVj) 



(<Pi\£ ai \<Pj) 



Xi % Xj 



Aa 2 

(-!)'+ ' 
Ao 2 



(1 - Xi)(l - Xj) 



iV(iV+ 1) - 



1 



A 



(1-Xi)(l-Xj) 



X iXj 



-N(N 



x 



j J 



(C15) 



(C16) 



(C17) 



(CIS) 



Although this does not appear clearly, one can verify that the matrix elements of 
operator T + £ a , 2 — £ ai are symmetric in accord with the fact that this operator is 
Hermitian over the region (01,02). 
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